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1. Introduction. It is the purpose of this paper to discuss the location of the roots 
of polynomials of high degree, with particular reference to the case of complex roots. 
This is a problem with which we at the Laboratories have been much concerned in 
recent years because of the fact that the problem arises rather frequently in the design 
of electrical networks. I shall not give any attention to strictly theoretical methods, 
such as the exact solution by elliptic or automorphic functions: nor to the develop- 
ment of roots in series or in continued fractions, though such methods exist and one 
at least—development of the coefficients of a quadratic factor—is of great value in 
improving the accuracy of roots once they are known with reasonable approximation. 

Instead, we shall deal with just two categories of solutions: one, the solution of 
the equations by a succession of rational operations, having for their purpose the 


dispersion of the roots; the other, a method depending on Cauchy’s theorem regarding 
the number of roots within a closed contour. 


PART I—MATRIX ITERATION 


2. Duncan and Collar. We shall treat the first category by a method recently 
elaborated by Duncan and Collar in two papers in the Philosophical Magazine.? I do 
not know how thoroughly these writers appreciate the close relationship of their work 
to that of the other writers whom I shall mention in the course of my presentation. 
The fact that their interest was primarily concerned with certain broad dynamical 
problems may perhaps have inhibited them from taking some of the steps which I 
shall take in their name. But they at least possessed the essential idea, and exhibited 
quite sufficient ability in the development of it to warrant the assertion that my 
presentation only differs from theirs in detail—sometimes details of omission, some- 
times details of amplification. 


* Received Dec. 26, 1944. 

1 The essence of this method is contained in a section of Legendre’s Essai sur la théorie des nombres. 
It is also attributed to Bairstow by Frazer and Duncan. It was developed independently, and perhaps 
somewhat more fully, by the present writer; but the extensions seem so obvious that it has not appeared 
to warrant separate publication. 

2 W. J. Duncan and A. R. Collar, A method for the solution of oscillation problems by matrices, Phil. 
Mag. (7) 17, 865-909 (1934); Matrices applied to the motions of damped systems, Phil. Mag. (7) 19, 197-219 


(1935). 
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3. The fundamental identity. We begin by noting that the A-determinant 





| @u +r da me's ae 
| @i2 Gog tA +++ Gne 
D(d) = }=[T@+ a) (1) 
| Gin den ci Gan + . 
is the characteristic function® of the matrix 
M = |\a;;|| (2) 
and its determinant 
A=| ai;|. (3) 


It is obviously a polynomial of degree n, which we may write 
D(A) = A" + pir™! — por”? +--+ + po. (4) 


Any quantity which satisfies the equation 


D(x) = 0 (5) 


and obeys the associative and commutative laws of algebra—whether it be a number 
or not—must also satisfy the relation 


A” = — pPr™! + poar™? +--+ F pa 
and if we multiply this by \ throughout and then eliminate \" we get 
n+1 2 n—1 n—2 
r = (pi + po)r — (pipo + ps)r ~ ote 5 


which is of the form 


n+1 (n+1)_ n—1 (n+1)_ n—2 
d = pi » + po r aM 
Similarly, by a continuation of the same process we may get a succession of equa- 
tions, all of the form 
m (m). n—1 (m)_ n—2 
hH =p dO +h2 dX +°°:. (6) 
We call the typical polynomial on the right of (6) f(A): graphically it represents a 
curve of degree n—1 passing through the m points —),;, (—A,)”. But we do not wish 
to emphasize this geometric interpretation but rather the formal algebraic fact that 
our derivation has required only the elementary rules of algebra and the relation (5), 


and that when these rules are satisfied 
A™ = f(A). (7) 
Suppose, now, that we expand the quotient f,(A)/D(A) in partial fractions. The 


result is 
fm(X) * fa(— dd 1 


Dis) X 4h; TH-~-k+e) 


j=1 
kj 


(8) 





3 The unconventional pecularities of sign in (1) and in (4) below happen to be convenient for our 





purposes later on. 
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But fm(—A;) =(—d,;)™ by (7), and DAA) = [J A+ A;) by (1). Hence 
Me+d 
ind) = Do (— 0) n(-— ). (9) 
kj At — dj 


Now (9) is an algebraic identity, and though we have used the process of division 
in setting it up, it does not require division by \ as a process of verification. Hence 
it is again true that if \ is any quantity which obeys the distributive and associative 
laws, such for example as a differential operator, and which satisfies (5), 


hm = D> (—d,)"w/ (A), (10) 
j=l 
where 
he 42X 
ri(d) = n(——). (11) 
hej NAR AG 


Note that the quantities denoted by 7;(A) are polynomials of degree m —1 in \ and are 
independent of m. 

4. Matrices. We next observe that, though matrix multiplication is not in general 
commutative, it is so if we restrict ourselves to certain groups. In particular, if we 
begin with the unit matrix J, any other matrix M, and all scalar quantities (i.e., 
numbers), then all matrices which can be formed from these by a finite number of 
additions or multiplications are commutative. For obviously M is commutative with 
itself and its powers, and with J, and with scalars, which observations together 
with the associative law are sufficient to warrant the general statement. 

Furthermore, we know from the Hamilton-Cayley Theorem‘ that 


D(M) = 


where D(A) represents, as in §3, the characteristic function of M. In other words, 
M satisfies all the requirements imposed upon \ in deriving the identity (10), whence 
we conclude that 


M™ = >> (—,)™*'(M), (12) 


where 7/(M) is a matrix independent of m. 
As a final step, we multiply this equation throughout by an arbitrary matrix 
K—which need not be commutative with the rest, since we shall perform no further 


operations—thus obtaining 
De (= Ay) 7 i(M) (13) 


where 7;(\/) = 7} (M)K is again independent of m. 

This is the fundamental identity upon which Duncan and Collar rely for their method. 
It is equivalent to n? equations of similar form connecting corresponding elements in 
the various matrices. For example, if a» is written for the element in the i-th row and 


M"K 


j-th column of MK, and e; for the correspondingly placed element in 7j(M/), it must 


be true that 


an = é(— di)” + €2(— Ao)” + lh sl oa én(— + aa (14) 


«M. Bocher, Introduction to higher algebra, MacMillan, New York, 1929, p. 296. 
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We again recall that —), is a root of the characteristic equation of —that is, 
of the polynomial D(\)—and hence is a number. The set —), are, in fact, just the 
roots which we wish to obtain. Similarly, the e;’s are numbers independent of m. But 


Q@m is a numerical function of m. 

5. The roots. Suppose now, that one of the roots which we will call —A,, is larger 
in absolute value than all the rest. Then if we select corresponding elements a,, and 
Qm+1 from two consecutive orders of MK we will have 


2 eo m+1 €3 Ag m+1 
1+— +—(—) +--- 
Am+1 4: \Ai é1 \Ay ‘ 


Om é2 f \2\™ és { \3\™ 
1+—(—) +—(—) +:-:: 
ej Ai é1 Ai 


and hence obviously 











Am 
ie ean ab he (15) 


mo Am 





In other words: if an arbitrary matrix K is multiplied repeatedly by M, and if its 
characteristic equation has a largest root, then the ratio of corresponding elements in two 
consecutive products approaches this largest root as a limit as m—>~. 

Similarly, we readily find that 


} Om+1 Om dA; m X; Xr; 2 
- ontare Bol) G/B, 0 
a Am—1 | , y™( : > : ArA2 dj ri 


whence if |Ai| and |)s| are greater than all other |)’s|, (whether they are them- 
selves equal or not), we again have 

















| Am+1 Om 
: | Am Am-1 | on 
lim = (— Ax)(— Ag). (17) 
m— © Am Am—1 | 
| 
Am—-1 A&m-2 | 
In the same way it can be shown§ that provided |Aa| poke my | A.| are all greater than 
[sss a [An 
| Am+i OQm+i-1° * *° @m+1 | 
Am+i-1 Am+i--2° °° Am | 
, Am-+1 Am = Qm—i+2 | 
lim — = (— Ai)(— Az) + + + (— Ad). (18) 
mM | @m+i—1 Am+i-2*** Am 
Am+i-2 Am+i-3** * Am—1 | 
| Qm Am—1 ge ti Am—i+1 





’ A. C. Aitken, Proc. Royal Soc. Edinburgh 46, 289-305 (1926), obtains formulae equivalent to 
these in a discussion of Bernoulli's method. 
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These equations are sufficient to determine all the roots in the particular case 
where 
[rr] >| Ao] >| As}--* >| An. 


6. Example. As a simple example we may take 


3 4 1 
wea ol ee loll: 
: a 0 
in which case 
3 11 43 
vea|2fee-|%f we-[ 
| 2 10 42 
171 683 
M‘K = I MK = | . 
170 682 


Taking the ratios of the first elements of consecutive matrices we get as the 
successive approximations to — [A] ; 


11/3 = 3.667, 43/11 = 3.909, 171/43 = 3.977, 683/171 = 3.992. 


Similarly we find that 








| 171 43 | | 683 171 | 
43 11 171 43 
= and = 4, 
| 43 11 | | 171 43 | 
i: .s 3s 


which should be the product of A; and dg. 
The characteristic equation in this case is, however, 


A+ 3 1 


D(a) = 
0) 2 A+ 2 


es 





and its roots are —1 and —4. The approximation is obvious. 

7. Complex roots. So far we have considered only real roots: for obviously, since 
complex roots occur in conjugate pairs (the coefficients being assumed to be real) 
there can be no Jargest one. Suppose, then, that |A:! =|\;| and that all other roots 
are smaller in absolute value. Then by (17), 


| 
Am+1 Am | 
| 
' 








— am QAm—1 





Qm-1 GAm-2 


This gives us the absolute value of the roots. It does not, however, determine the 
angles. To get this, we can best return to equation (14) and write (retaining only the 


leading terms) 
Qm = €1(— A1)™ + €2(— Az)”. 
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Writing the similar equations for m—1 and m-+1, and eliminating e; and @, we get 
A1A2@m-1 + (Ay T No) Qm + Am+1 = 0. 
Substituting the value of \iA2 as given by (19) we get finally 
Am +1 Am—1 
. | Qm Am—2 
—-i+') = in —————. (20) 
m—> @ Am Qm—1 | 
Qm—1 Am—2 
This, together with (19) is sufficient to determine the pair of roots. 
As written, the formula applies even if the roots are real.6 When they are complex 
it is best to write —A,;= —\_,=pe‘*. Then obviously we need only replace the dj): of 


(19) by p?, and the —(Ai;-+Az) of (20) by 29 cos ¢. 


Similar, but more complicated, formulae can be obtained when more than two 


roots have the same absolute value. 
8. The method of Daniel Bernoulli. We now note that any polynomial in A, which 


we take in the form 


D(d) = A* + pit”! — pod"? +--+ F paid + Pn (4) 


as before, can be written as 


| A 1 0 0 0 
}O nr 1 0 0 | 
-_ ££ * 0 0 | 
D(dr) = | oh | (21) 
0 0 0 r 1 
Pn = Pn—-1 Pn—2 pe pith 
But this is the characteristic function of the matrix 
| 0 1 0 sore K@ 0 | 
| °° 0 tf +. @ @ 
ae as } 0 0 0 «+. O 0 | (2) 
| 0 0 GS «. 2 1 
Pn pn-1 Pao'** po fpr 


Hence if we choose for K any matrix whatever, we may solve for the largest roots by 
any of the equations of §§4 and 6. 

It is particularly convenient to take K in the form 

* It is not even necessary that they be equal in absolute value, though unless they are equal (or 


nearly equal) (15) will obviously be a more convenient formula. 
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10 0 ao 
0 0 ai 
K=1|10 0 a2 (23) 











0 0 O-+* ani 


Then we have 











0 0 O-++ ay, 
0 0 0O-- ae 
0 0 O-:- a; 
MK = ’ 
n—l1 
| 0 0 O--+ Do pa jors| 
j=0 





which is again of the same form as K. If we denote ) 20 pj+10n—j-1 by an, we also have 


M*K =||}0 0 O:-:-a@% 


n—l 
Oo 8 6. « > Pj-10n—j 


j=0 














And in general 
° An 


* Am+1 


M"K =||0 0 0--- amis 














0 0 O¢+ + usa 
where 
n—1 
ak = ye P j+10K— jy k>n-— 2. (24) 
j=0 
This entire set of matrices, however, is characterized by a simple sequence of a’s, 
of which the defining equation is (24). Obviously, it is also true that any set of four 
consecutive a's in this sequence also constitutes a set of corresponding elements from four 
consecutive matrices of the set M"K. Hence, the use of the symbol a@ in this connection 
is consistent with its use in §§3—6. But (24) is the recursion formula used in Ber- 
noulli’s method of solution as developed by Euler, Lagrange and Aitken. Hence this 
particular special case of the results of Duncan and Collar is identical with Ber- 
noulli’s method. 
Concerning this method Whittaker and Robinson? say: “Though hardly now of 
first-rate importance, it is interesting and worthy of mention.” Our tests at the 


7 E. T. Whittaker and G. Robinson, Calculus of observations, 2nd ed., Blackie & Son, London, 1929. 
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Laboratories, however, have shown it as good as any other method in the case of 
complex roots. Such inferiority as it may have compared to the root-squaring method 
as regards speed is quite compensated by the fact that it is self-correcting: that is, 
an error at any stage of the process merely prolongs the calculations, but does not 
invalidate it. 

9. The method of R. L. Dietzold. Another form into which the general results of 
Duncan and Collar can be thrown is obtained by using the conjugate form of (21) 
together with the same matrix for K as before. Denoting the conjugate of M by M’, 
we have from (22) and (23) 


eee T 
0 0 O oe ao + Pn-1%n-1 
M’'K = || 0 0 0 ae 2 Qa) + Pn—24n-1 | 
| | 
| 0 0 0 a An—2 a P1Qn-1 | 
If, then, we define 
ag = PnOn-1, a; = Aj-1 + Pn—jOn-1; (25) 


M’'K becomes identical with (23), except that all the a’s are primed. In general, if 


we set 
(m) (m—1) (m—1) 
(26) 


Qj = api t+ Pn—jOn-1 ; 
and understand that a”? is zero for all m, we have 


(m) 
0 0 0 2 * 
0 


(m) 


| 
| 0 O--++a, 1] 
u"K =| 6 O:-+e,” Ih (27) 
| | 
| | 
(m) 


10 0 O+ ++ aa_-1 |} 


In this case, as in all others, the index m is the one which is to be varied in using 
formulae such as (16)—(20). 

This variant of the general scheme of Duncan and Collar was developed by Mr. 
R. L. Dietzold of the Bell Telephone Laboratories, but has not been published. 
As compared with Bernoulli's, it has the merit of using a large number of simple 
operations instead of a small number of complicated ones. It is approximately as fast, 
and like all schemes based on Duncan and Collar’s results, it is self-correcting. 

10. The method of Graeffe. There is also a close connection between Duncan and 
Collar’s processes and the root-squaring method. This method, which is usually 
attributed to Graeffe, seems actually to have been developed first by Dandelin, and 
has had the attention of a long list of mathematicians, including Lobachevski, Encke, 


Brodetsky and Smead, and Hutchinson. 
This connection can best be established® by recalling that the roots —A; of the 


8 M. Bocher, Introduction to higher algebra, MacMillan, New York, 1929, p. 283, Theorem 3. 
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matrix AJ+M are invariant under transformations of the type 7—[AI+M]T. 
Furthermore, it is possible to find a transformation of this sort which will throw M 


into the form 
+>. On 


| 
Se hk-+-4 
M* = TMT = 
6 On-irk& 
and hence \J+ M into the form A\J+M*, since T-YT is obviously J. 
This same transformation, however, carries M™ into M*", as we readily see from 


the identity 


T-\(MTT-MT - ++ TM)T 
(T7MT)(T-MT) - ++ (TMT) 
= M*., 


| ie ies i 


Ih 


Hence the characteristic equations of M™ and M*” must also have identical roots. 
But, obviously, 

m 

1 


m 


|| » 0 
aa er 
| 





M 


ie 8 -*2 kh, 


so that the roots of its characteristic equation, and therefore also those of the char- 
acteristic equation of M™, must be —Aj’. 

But if we take K=J in Duncan and Collar’s process of matrix iteration, the suc- 
cessive matrices obtained are M™. Hence the whole process may be regarded as one 
which sets up a sequence of characteristic equations with roots —\;, —Aj, -- - and 


in general —)j’. 
In the root-squaring process as originally developed only the powers —\j, —Aj, 
—*,--- were obtained, which corresponds in matrix terms to getting first the 


product of M by M, which is M?; then the product of M? by M? which is M‘, and so 
on. Thus high powers are reached with a smaller number of matrix operations, which 
is theoretically desirable. Practically, however, the superiority is not so apparent. 
For the zeros of (22) are rapidly replaced by numbers in forming powers of M, so 
that a multiplication such as M*-M® involves many more arithmetical operations 
than a multiplication of the form M-M®*. Furthermore, an error at any point of the 
root-squaring method perpetuates itself, whereas in the other method an error at 
any stage is merely equivalent to starting over again with a new value of K. 

Our experience leads us to believe that the methods of §§8 and 9 are generally 
to be preferred, at least when computations are to be performed by a clerical staff of 
computers. 

11. The method of Bernoulli as developed by Lagrange. There is also a very close 
connection between the iterated matrix M™ and a development of Lagrange’s which 
he characterizes as based upon that of Daniel Bernoulli. In it, he notes that 
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D’(d) 1 n Sy Se S3 
D0) ity tw tT BRT (28) 
where 
Su (= Aa)* + (— Ae)” + - >> + (— A)"; (29) 


and it is the quotient s»/sSm—1 which Lagrange uses. Obviously, these are just the sums 
of the elements in the principal diagonals of M*”. But Lagrange’s method of obtaining 
them by dividing D(A) into its derivative is preferable. Besides, in spite of what 
might at first be assumed, it is self-correcting. 

It is of historical interest to note that a very similar development was worked out 
by Legendre® independently of Lagrange, and at about the same time. Both of these 
writers, however, knew of earlier work by Euler, who had carried out a similar de- 
velopment using instead of D’(A) an arbitrary polynomial of degree n—1, which 


A- PLANE D- PLANE 


c c 


, ) 
a Od 














Fic. 1. 


merely has the effect of replacing the s,’s in the right-hand member of (28) by the 
Qm's defined by (14). In other words, the method of Euler was exactly equivalent to 
that of Duncan and Collar, except that in the former there was no obvious criterion 
for the choice of a convenient form of numerator, whereas it is easy to choose matrices 
K which will lead to a simple succession of operations, as we have illustrated in 
Sections 8 and 9. 

PART II—CONFORMAL MAPPING 


12. The method of Routh. The second group of methods to which I wish to refer 
are all founded upon a well-known theorem of Cauchy. If we represent the complex 
variable \ by one plane, and the complex variable D by another, then the equation 


D(X) = A* + prr*"! — por”? +--+ + pny (4) 


may be looked upon as a transformation by means of which the A-plane is mapped 
upon the D-plane. The correspondence between \ and D, however, is not 1:1 but in 
general 2:1; and hence a simple closed curve C in the A-plane (Fig. 1) passes into a 
much more complicated curve C’ in the D-plane. In regard to the curve C’ the the- 
orem in question says that the number of times it loops around the origin is exactly 
equal to the number of roots of D(A) =0 which lie inside C. 





9 Legendre’s development was in terms of the reciprocal powers of the roots, instead of their direct 


powers. Otherwise the two were identical. 
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This rule appears first to have been applied by Routh to the problem of determin- 
ing the number of roots with positive real parts, a problem which interested him be- 
cause of its relation to the stability of linear dynamical systems. For this purpose he 
used as the contour in the \-plane the imaginary axis closed by a semicircle of infinite 
radius, thus enclosing the entire right half of the plane. For this particular contour 
he explained in great detail how from the sequence of intersections of C’ with the 
real and imaginary axes the number of roots could be found without more definite 
information as to the shape of C’. He also developed a sequence of functions, similar 
to Sturm functions, by means of which the number of roots could be determined 
from the polynomial directly without even knowing the real and imaginary inter- 
cepts of C’. He did not extend either of these studies to the point of locating the roots 
more exactly, but both are capable of such extension and have actually been used. 

13. The method of G. R. Stibitz. The second method—the one using functions 
similar to the Sturm functions—was developed further by G. R. Stibitz of the Bell 
Telephone Laboratories. He observes, first, that the method can also be used to find 
the number of roots with real parts greater than Xo. To do this, it is merely necessary 
to replace \ by A—Ao in the polynomial (4), and then proceed as outlined by Routh. 
By carrying out this process for enough values of Xo, the roots can be segregated within 
strips parallel to the imaginary axis. Then by a definite routine (resembling in its 
essentials the Weierstrass subdivision process in point-set theory) the real values of 
the roots can be found to any desired degree of approximation. When this has been 
accomplished, the imaginary parts are determined at once as a ratio of two of the 
Sturm-like functions. 

Stibitz has developed complete schedules for the computations required in solving 
polynomials by this method, for all values of m up to 10. The method has been tried, 
and works reasonably well, though perhaps not as rapidly as those explained in 
Sections 8 and 9. I suspect that the decision in this case, however, must remain a 
conditional one; for the computational routine of Stibitz’ method is complicated 
(i.e., varied) as compared with the extremely simple (i.e., repetitive) routines of 
Sections 8 and 9. For this reason, it is not as well adapted to use in an industrial 
computing laboratory. In the hands of a mathematician who thoroughly understood 
its theoretical origin it might show up much better. 

14. The method of A. J. Kempner. Kempner’s methods" resemble more nearly the 
other portion of Routh’s work. He chooses as his contour C a circle of radius r about 
the origin as center. Then \ = re, and (4) becomes 


D(d) = [r" cos n6 + pir"! cos (n — 1)0 — por"? cos (n — 2)0+--- | 


(30) 
+ ilr™ sin nO + pyr™ sin (n — 1)0 — por™? sin (n — 2)0+--- |. 


Thus the real and imaginary parts of D are trigonometric sums, which, as Kempner 
remarks, could be calculated by means of a harmonic synthesizer, such for example 
as the Michelson “analyzer.” Thus two curves would be obtained, one giving the 
real part of D, and the other its imaginary part, both as functions of 8. From this 
point on, Kempner suggests two possible routines. First, to regard these curves as 
parametric representations of D, and from them construct the curve C’ itself. Second, 
to keep them as separate curves and not bother further about C’. In both cases, he 


1° University of Colorado Studies 16, 75 (1928); Bulletin of the Amer. Math. Society 41, 809 (1935). 
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develops rules very similar to Routh’s for finding the number of roots directly from 
the sequence of intercepts with the axes. 

He uses this routine to segregate the roots in annular rings, and then tracks down 
their absolute values by a suitably chosen succession of intermediate circles. The 
angle of any root is, of course, automatically determined as the value of @ at which 
the real and imaginary parts of (30) vanish when r is given the particular value 
appropriate to that root. 

Kempner also mentions the possibility of applying the method to sectorial in- 
stead of annular regions, but does not develop this idea to a significant degree. 























SASSY 

















Fic. 2. 


15. The isograph. Kempner’s method was also developed independently, but 
somewhat later (1934) at Bell Telephone Laboratories, and led to the construction 
of a machine, called the isograph, which draws the curve C’ corresponding to a circle 
of any radius r. 

Since the independent variable in plotting the curves is an angle, what is required 
for the isograph is a rotating unit that provides two linear motions—one proportional 
to the sine and the other to the cosine of the angle. There would have to be ten of 
these units to provide for the ten variable terms of a tenth degree equation, and while 
the first unit moves through an angle 6, the second unit must move through an angle 
20, the third unit through an angle 30, and so on. Then by providing a means of 
summing the sine and cosine motions separately, and allowing these sums to control 
two perpendicular motions of a pencil and drawing board, a closed curve will be de- 
scribed as @ increased from 0 to 360 degrees. 

To secure motions proportional to the sine and cosine of the angle of rotation, 
the isograph utilizes the “pin and slot” mechanism illustrated in Fig. 2. Here an arm 
rotating about a fixed point carries a pin arranged to slide, by means of a rectangular 
block, in rectangular slots cut in two slide-bars, each of which is free to move back 
and forth in one direction only—the two motions being at right angles to each other. 
These motions are equal to the length R of the arm times the sine and cosine of the 
angle of rotation. 
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The ten units provided are geared to a common driving motor, but the gearing is 
designed so that when the arm of the first unit moves through an angle 0, that of the 
second unit will move through an angle 26, that of the third through 30, and so on. 

To provide for summing up all the sine terms and all the cosine terms, the ends 
of all the slide-bars carry pulleys so that a single wire may be carried around all the 
sine pulleys and another around all the cosine pulleys as indicated in Fig. 3. Station- 
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ary pulleys are mounted between the movable ones so as to keep the direction of 
pull on the wires in line with the motion of the slide-bars. These wires control the 
relative motions of a pencil and drawing board to plot a curve as the angle is varied 
from zero to three hundred and sixty degrees. 

The construction of the rotating elements is shown in Fig. 4. The drive shaft 
passes through the bed plate and is fastened to the center of a steel bar that acts as 
the arm of Fig. 2. This bar is grooved to receive the pin of the “pin and slot” mecha- 
nism. In order that the pin may be adjusted for different crank lengths, corresponding 
to the coefficients py7,-, of the various terms in the equation, a rack is cut along 
one edge of the groove so that a pinion attached to the pin may move it along the bar. 
After adjustment the pin is secured in place by a set-screw. 

The top of the bar carries a carefully graduated scale to which the center of the 
pin must be set accurately. The scale is made visible at the center of the pin by con- 
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structing the latter as a hollow cylinder. A vernier scale within the cylinder enables 
the effective arm length to be adjusted very exactly to the desired value on either 
side of the center—one side for positive coefficients and the other for negative. The 
total range of adjustment is three inches. 

The hollow pin turns in a rectangular bronze block which fits the slots of two 
slide bars, one for the sine motion and one for the cosine motion. The slide bars are 
identical steel plates running in bronze ways set accurately at right angles to each 








Fic. 4. 


other. At the end opposite to the slot each plate carries a pulley around which is 
passed the wire that sums up the sine or cosine motions of the ten elements. One 
end of each wire is fixed. The other end of the cosine wire is led by pulleys to the 
drawing board, which consists of a thin aluminum sheet mounted on ball-bearing 
rollers so that it is free to move back and forth in only one direction. A counterweight 
fastened to the other edge of the board keeps the wire under constant tension. The 
free end of the sine wire is led by pulleys to a counterweighted pencil carriage, which 
is mounted with ball bearings in a fixed guide crossing the drawing board at right 
angles to its direction of motion. Thus the board is displaced back and forth in 
proportion to the sum of the cosine terms, and the pencil is displaced back and forth 
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in a perpendicular direction in proportion to the sum of the sine terms; and this com- 
bined motion gives the desired curve. 

In operation, the isograph has given accuracies of one per cent or better; and of 
course gives them quite rapidly. In fact, the most rapid method we have at present 
is that of using the isograph to obtain this degree of accuracy, and then improving 
it either by the methods explained in §§8 and 9, or by successive approximation to 
the quadratic factors. 

16. Conclusion. In conclusion I wish merely to point out that in none of the 
methods which I have described is computation with complex numbers involved. 
They are all real methods. At present this seems to be a fundamental requirement 
imposed upon us by commercial computing machines, since the multiplication of two 
complex numbers on such machines requires six, and division eight, separate opera- 
tions. If this restriction were removed, other methods might conceivably prove to be 
more rapid. 

Partly with this in mind, and partly because we must frequently deal with com- 
plex quantities in other connections, we are at present developing a computing 
machine for complex quantities. When it is completed, as we hope it will be in the 
course of the present year, we shall undertake a further study of methods which now 
are clearly ruled out by mechanical limitations.* 

POSTSCRIPT BY R. L. DIETZOLD 

When the foregoing paper was written, it was intended for immediate publication. 
By coincidence, however, several other papers of similar character appeared at just 
about that time, and Dr. Fry concluded that the subject was of too limited interest 
to justify publishing another. 

Since then, the situation has changed in several ways. First, the interest in meth- 
ods of numerical computation has greatly increased, largely because war activities 
have led to much work of that kind. Second, the specific problem of root-finding has 
become a live one because its fundamental importance in linear dynamics is more 
widely recognized. Finally, a few new methods of iteration have been evolved and 
some new types of computing machines developed. The paper therefore now has a 
timeliness which it lacked when written, but a few comments are required to bring 
it up to date. The most important of these are noted in the following paragraphs. 

In the Bell Telephone Laboratories the available computing equipment has been 
materially improved through the development of the relay computer by Stibitz and 
this inevitably reacts upon the relative convenience of various methods of solution. 
Although the relay computer is very flexible in respect to the type of problem it can 
handle, it is particularly well suited to iterative processes such as Bernoulli’s method 
of root extraction; for once the proper instructions have been set into the control 
tape which governs the machine, all successive operations are performed without 
further supervision. The simplicity of Bernoulli's rule, which requires only that the 
machine accumulate »—1 of the a’s, each multiplied by the appropriate coefficient 
from the polynomial, recommends it for mechanization. The instructions are easily 





* This machine was placed in service in 1940 and was demonstrated at the summer meeting of the 
American Mathematical Society in Hanover, New Hampshire in September of that year. The relay com- 
puters referred to in Dr. Dietzold's postscript are still more versatile devices which have been developed 
since that time. 
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set up, and the machine is not required to recall very many numbers at any stage in 
the process. Bernoulli’s method is likewise well adapted to computing equipment of 
the punched-card type, provided only that the accumulator is designed to recognize 
algebraic sign. 

One of the routines which may be set up in a relay computer enables the algebraic 
operations to be performed on complex numbers with the ease that the same opera- 
tions are performed on real numbers with a mechanical computing machine. The 
availability of this aid makes Newton’s method useful for root improvement in the 
complex domain, and some on the Laboratories’ computing staff prefer it to Bair- 
stow’s method, although the margin of choice is not great. 

Bairstow’s variation of Newton’s method avoids computation with complex 
quantities by improving the coefficients of a trial quadratic factor. The trial factor, 
say Q(A) =\?+a\+), is divided twice into the polynomial, and the rates of change of 
the remainder coefficients found from the second remainder, as in Horner's process. 
The method has by now been sufficiently publicized ;“ nevertheless, it can be given 
here, since it is short to state. The polynomial being expressed as 


D(d) = (rod + 50) + Q(A)(riA + 51) + O(A)(to + AA+---), 


improved coefficients for Q are 


| To ri | ar; — S$, fo | 
: So $1 | br; So 
¢ = 36ers" : b + ——_—__ 


| Oy — $3 r; | | Gy = 33 r, | 





| 
7 


| bry $1 | br; 51 | 

Newton’s method typifies a class which is deliberately excepted from treatment 
in Fry’s paper; methods in this class are characterized by the property that only some- 
times do they lead to a solution. Newton’s method, for example, can never lead toa 
complex root if the iterative process is started from a real trial value. Bairstow’s 
method has a similarly restricted region of convergence and was, quite properly, 
advanced by him only as a means for improving roots already located approximately. 

Methods which sometimes fail to converge may still be very useful if, in appli- 
cation, they converge often enough and fast enough. Newton’s method and its varia- 
tions, however, almost always fail unless they can be started from values closely 
corresponding to roots. But in 1941, Shih-Nge Lin revealed an algorithm” remarkable 


1! Bairstow gave the method only in Reports and Memoranda No. 154, Advisory Committee for 
Aeronautics, Oct., 1914 (H. M. Stationer’s Office), but it was made generally available by Frazer and Dun- 
can, Proc. Royal Soc. London 125, 68-82 (1929). Hitchcock offered the method as An improvement on 
the G.C.D. method for complex roots, Jour. Math. Phys. 23, 69-74 (1944). Hitchcock proposes that the 
roots be improved by this method after only approximate location by the G.C.D. method, which he gave 
in Jour. Math. Phys. 17, 55-58 (1938). The G.C.D. method is nearly identical with the method of G. R. 
Stibitz, described by Fry. Bairstow’s method was also rediscovered by Friedman, whose work is noted in 
Bull. Amer. Math. Soc. 49, 859-860 (1943). Bairstow’s formulae give the leading terms of series develop- 
ments of the coefficients by Fry, who concluded, after an investigation of the convergence, that the ex- 


pansion was suitable only for root-improvement. 
2 A method of successive approximations of evaluating’the real and complex roots of cubic and higher 


order equations, Jour. Math. Phys. 20, 153 (1941). 
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both for simplicity and canvergence. By’Lin’s method, the polynomial is divided 
only once by a trial quadratic factor; if® 

pot prt pr+--- 

(rod + So) + O(A) (go + gid + god” +- A F 


D(\) 


Il 


II 


improved coefficients for Q are 
; gop: — Qipo , Po 
a’ = —————_ i =—. 
1% 19 
If it converges, the process determines the factor corresponding to the roots of 
least absolute value; thus a suitable initial choice for Q is 


dX? + (pi/ po) + (po/ po). 


In application, the process does very often converge, although sometimes slowly. 
When the convergence of Lin’s method is slow, Bairstow’s method offers a valuable 
supplement. Lin's method is used until the size of the remainder indicates that an 
approximation to a quadratic factor has been obtained; Bairstow’s process, started 
from a sufficiently close approximation, will converge, and when it converges, it 
converges rapidly. . 

The combination of these two methods provides useful, and usually adequate, 
equipment for the work-a-day solution of polynomial equations. In recalcitrant cases, 
mechanical aids are particularly helpful. Bernoulli’s method is always available, but 
is quite likely to be slow in cases for which Lin's method has already failed. This 
makes little difference if the iterative process is performed automatically by a relay 
computer, but recommends devices to accelerate the convergence if the computation 
must be performed without aid. An efficient device for accomplishing this is given by 
A. C. Aitken in a very full discussion“ of numerical methods for evaluating the 
latent roots of matrices. 

Like most of those who use matrix methods, Aitken is concerned not solely with 
the solution of polynomial equations, but rather with the more general problem of 
determining the characteristic roots (and also the characteristic vectors) of matrices. 
Preliminary reduction of the matrix to the rational canonical form involves so many 
operations,” that one would commonly start the general problem with a matrix M 
having few vanishing elements. In this event, we lose one of the reasons for preferring 
Bernoulli's method (i.e., repeated multiplication by M) to matrix powering by the 
root-squaring method, for the latter method arrives at high powers of M with fewer 
operations, thus providing another means for hastening the convergence. The ad- 
vantage is, however, partly illusory except for the limited class of computers who are 
so unerring that they can afford to sacrifice the self-correcting feature of the former 


procedure. 


43 A departure from Fry’s notation is convenient here. 
™ Proc. Royal Soc. Edinburgh 57, 172-181 (1937). 
5 Harold Wayland, Expansion of determinantal equations into polynomial form, Quarterly Appl. 


Math. 2, 277-306 (1945). 
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THE KARMAN-TSIEN PRESSURE-VOLUME RELATION IN THE 
TWO-DIMENSIONAL SUPERSONIC FLOW OF 
COMPRESSIBLE FLUIDS* 


BY 


N. COBURN 
University of Texas 


1, Introduction. T..v. Karman and H. S. Tsien! have treated the two-dimensional 
subsonic flow of a perfect, irrotational, compressible fluid by replacing the adiabatic 
pressure-volume curve by the tangent line drawn at an arbitrary point of this curve. 

First, we shall discuss the applicability of the Karm4n-Tsien idea in the supersonic 
range. Secondly, we shall show that when the Kérmén-Tsien relation can be used 
(fairly uniform completely supersonic flow), the characteristics form a Tschebyscheff 
net (fish net).2 However, we shall be concerned with those regions of the physical 
plane which can be mapped into a Tschebyscheff net in a unique one-to-one manner. 
Hence, we shall not study the onset of shock. Further, we shall show that if the di- 
agonal curves of the net of characteristics are drawn so as to correspond to equi- 
distant values of the arc length parameter along the characteristics, then these 
diagonal curves will be the families of equipotentials and stream lines. Analytically, 
this last result means that the determination of the stream lines depends upon two 
arbitrary functions of one real variable. It is shown that the angle between the char- 
acteristics and the angle formed by a tangent to a stream line and the x-axis can be 
determined in terms of these functions. Further, the magnitude of the velocity and 
the density depend upon only the former angle and the Mach number of the flow. 
In particular, if a known stream line coincides with the x-axis, it is shown that only 
one arbitrary function enters into the problem of determining the stream lines. Even 
in this last case where the data are of a simple Dirichlet type (symmetric flow about 
the x-axis and a known external boundary stream line—as in the jet problem), the 
direct problem cannot be solved easily. Hence, an analytical-geometrical method is 
outlined for solving the problem indirectly. A particular example is studied. Finally, 
in an appendix, we furnish another proof (analytical) of the fact that when the 
K4rm4n-Tsien relation is applicable, the characteristics form a Tschebyscheff net and 





conversely. 
2. Extension of the Kaérmén-Tsien method to supersonic flow. In this section, we 


shall show that the Kaérman-Tsien method may be extended to the supersonic flow of 
a perfect, irrotational, compressible fluid. If we denote the pressure by p, the density 
by p, the ratio of the specific heats by y, the adiabatic relation is 


pp-? = constant. (2.1) 


* Received Oct. 16, 1944. 
1T. von Karman, Compressibility effects in aerodynamics, Journal of Aeron. Sciences 8, 337-356 


(1941). 
H. S. Tsien, Two-dimensional subsonic flows of compressible fluids, Journal of Aeron. Sciences 6, 


399-407 (1939). 
2 L. Bianchi, Lezioni di geometria differentiale, vol. 1, Enrico Spoerri, Pisa, 1922, pp. 153-162. 
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Replacing the isentropic curve (2.1) by the tangent line drawn at the point (p7", pi) 
in the pressure-volume diagram (or by a hyperbola drawn at the point (1, p:) in the 
pressure-density diagram), Tsien obtains the relation 


bi—p=apilp —p), (2.2) 


where aj is the velocity of sound corresponding to (91, ~:). By use of (2.2), it is easily 
shown that the Bernoulli relation becomes (where w is the velocity) 


$s, - a 
o = wi = aipi(p F = .9 (2.3) 
Further, by use of (2.2) and the definition of a* (that is, a?=dp/dp), it follows that 


op @émal, (2.4) 


where & is some constant. Hence, (2.3) can be transformed into the following forms: 


y\2 ane 2 
GS -(*)- 1, (2.5) w—a=w—a=l, (2.6) 
ay a Pp 


where / is some constant. 

In the following, we shall assume that the point (p:, p:1) corresponds to a super- 
sonic state of the fluid. From (2.5), we see that as w increases, p decreases. Further, 
from (2.4), we see that as p decreases, a increases. As noted by Tsien, the first result 
is in accord with the physical facts; the second result is in contradiction to known 
physical facts. However, (2.6) furnishes some useful information. Since the density p; 
corresponds to a supersonic state of the fluid, the equation (2.1) is valid for this p; 
and the corresponding ;. Hence, by well known results, w,; is larger than a;. Thus 
from (2.6), we see that w is always larger than a. That is, the fluid is always in a super- 
sonic state in this sense of the term. However, by dividing (2.6) by a? and noting that 
as w increases, p decreases, and a increases, we see that as w increases, w/a decreases. 
This ratio approaches the limiting value 1 as w tends to infinity. Hence, the behavior 
of w/a is contrary to that of a real fluid. 

Perhaps the best indication of the permissible values of w which can be used for 
a given w; is obtained by following the procedure of Tsien. If we consider the upper 
limit of the useful values of w as occurring for p=0, we find that the corresponding p 


is given by 





2 
1 
1 ft oi (2.7) 


2,2 
p a*p' 


Substituting this value of p into (2.3), we obtain the relation 


(=)- fb eae I( = 1) - | (2.8) 
w,) (wi/ai*) L\ atp, 


Since the values of p1, p; satisfy (2.1), we find 


d 
dp/ Pl 
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w\? 1 1 " 
(=) =14 -- {(—+1)-1]. (2.10) 
Wy (w/a)? ¥ 


If y is taken as 1.4, a simple computation reveals that as w,/a; goes from 1 to ~, the 
ratio w/w, runs from 1.7 to 1. That is, for large values of the ratio w,/a,, the range of 
applicability of formula (2.2) is severely restricted as regards the upper limit of w. 
Hence, the Karm4n-Tsien relation should be useful in the supersonic range for a 
fairly uniform fluid flow. Further, as we shall show in the next section, the character- 
istics in this case form a Tschebyscheff (fish) net. We shall not be concerned with the 


Thus (2.8) becomes 


onset of shocks. 

3. The geometry of the characteristics for the relation (2.2). If u(x, y) and v(x, y) 
denote, respectively, the x and y components of the velocity for the steady flow of a 
fluid at any point P of the plane region considered and p(x, y) denotes the density of 
the fluid at P, then from the equation of continuity it follows that a stream function 
(x, y) exists such that 

ay ay 
(3.1) 


pu = —-> pu = —- 
dy Ox 


Further, since the motion is irrotational, a velocity potential exists such that 


“ 0d 0 (3.2) 


ax dy 
For a given pressure-density relation, the Bernoulli relation determines p as a func- 
tion of u?+v?. Hence (3.1), (3.2) constitute a non-linear system. Eliminating the par- 
tial derivatives of p from the continuity relation by use of the Euler equations, 
we find that $(x, y) satisfies* 
0° 0° 0° 


(a? — 4*) — — 2u0 —— + (a* — 2°) 


— = 0. (3.3) 
Ox? Oxdy dy? 


In the supersonic range, the equation (3.3) is hyperbolic. Let us denote the equa- 
tions of the two parameter family of characteristics of (3.3) by 


x = x(a, B), y = y(a, B), (3.4) 


where a=constant and 8 =constant are the parametric equations of the characteris- 
tics. We denote the arc length element of the net formed by the characteristics by 


ds? = A*(da)? + B*(dB)? + Cdadf, (3.5) 


where A?2, B?, C are the metric coefficients of the net. It follows from (3.3) that the 
projections of the velocity vector on the normals to the characteristics have the 
magnitude a. This means that the projections of a vector, normal to the velocity 
vector and of magnitude equal to that of the velocity vector, on the tangents to the 
characteristics have the magnitude a. Also, the projections of the velocity vector on 
the tangents to the characteristics have the magnitude \/u?+v?—a?. Further, from 


2? R. von Mises and K. O. Friedrichs, Fluid dynamics, Brown University, Providence, R. I., 1941, 


p. 230. 
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(3.1), it follows that p~ times the gradient of is a vector, normal to the velocity 
vector and of magnitude equal to that of the velocity vector; and, from (3.2), it fol- 
lows that the gradient of ¢ is the velocity vector. From the above properties of the 
characteristics and those of the gradient, it follows that 





oy ay 
— = pad, — = — paB, 3.6 
re 36 p (3.6) 
Od eee eae Oo in ae 

—=1/uv2+ 0? -— aA, — = Vu? + v? — a? B. (3.7) 
0a 0p 


We shall prove that the net of characteristics forms a Tschebyscheff net, when the 
Kdrmén-Tsien relation is applicable. 
From (2.4), (2.6) and (3.6), (3.7), we find 


Oy } 0 0 
ra, va a, (3.8) ehh, =k (3.9) 


Ja 0B da 0p 


The integrability conditions for (3.8), (3.9) furnish the result 


OA OB 
aoont teh eo tet (3.10) 

op Oa 
Hence A and B are functions of a and 8, respectively. By proper choice of scale fac- 
tors, A and B may be assigned the value unity. The new parameters a and £ are then 
arc length parameters and the net is a Tschebyscheff net. 

Next, we shall derive a result similar to that obtained by von Mises‘ in plane 
plasticity: when the Karmén-Tsien relation is applicable and the diagonal curves of the 
characteristics are drawn so as to correspond to equi-distant values of the arc length pa- 
rameter along the characteristics, then these diagonal curves will be the families of equt- 
potentials and stream lines. 

Since (2.2) is valid, we see from our previous result that a and 6 may be chosen 
as arc length parameters. Hence, the nets 

a+ 6 = 2 = constant, a — B = 2n = constant, (3.11) 
represent, respectively, the diagonal curves of the net of characteristics, correspond- 
ing to equi-intervaled values of the arc length parameters a and 8. Further, from (3.8), 
(3.9), we obtain 

icles sali oe (3.12) 
l k 
Hence, the diagonal curves =constant and y=constant represent, respectively, the 
equipotentials and stream lines. 

With the aid of our previous results and known properties of Tschebyscheff nets, 
we obtain some additional results. The general representation of the stream lines in the 
supersonic range for the Kaérmdn-Tsien relation depends upon two real arbitrary func- 
tions. If one stream line coincides with the x-axis, these two functions are equal except for 


‘R. von Mises, Bemerkung zur Formulierung des mathematischen Problems der Plastizitatstheorie, 
Zeitschr. fiir angew. Math. u. Mech., 5, 147-149 (1925). 
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a constant. Further, the velocity and density depend only upon the angle between the char- 


acteristics and the Mach number of the flow. 
For a Tschebyscheff net, it is well known? that (3.5) may be written in the form 


ds? = (da)? + 2 cos wdadB + (df)?, (3.13) 


where w is the angle between the two families of characteristics of the net at any point. 
Further, it is known that w may be expressed in terms of two arbitrary functions F(a) 
and G(8) by the relation 


w = F(a) + G(8). (3.14) 
Finally, the general representation of the net is given by 
a B 
t= f cos F(é)dt + f cos G(é)dt, (3.15) 
a B 
y= f sin F(t)dt -f sin G(é)dt. (3.16) 


Introducing the parameters along the equipotentials and stream lines from (3.11), 
we find that the above equations become 


w w 

ds? = 4 cos? ry (dé)? + 4 sin? ry (dn)?, (3.17) 

w = Fé +n) +G(E— 4), (3.18) 
E+ §—1 

z= f cos F(t)dt +f cos G(é)dt, (3.19) 
+9 &—-7 

y= f sin F(t)dt — f sin G(é)dt. (3.20) 


Another relation of the form (3.18) can be obtained by introducing the angle 
6(~, ») which the tangents to the stream lines (7=constant) form with the x-axis. 
Let the equations of a stream line be 


x = 2(s), y = y(s), n = constant, (3.21) 


where s is the arc length parameter along the stream line. From (3.19), we find by 
differentiation 
dx dé 
— = [cos F(é + ) + cos G(é — n)] ty (3.22) 
ds 


ds 
By use of the well known addition formulas of trigonometry, (3.22) becomes 


dx ss eas [= +n) +G(é —- 4 ae [= +) -Gg- a (3.23) 
ds 2 2 ds 








From (3.17), we find that along 7 =constant 


de 1 ‘ [a + GE — 4 
- 


(3.24) 
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Substituting (3.24) into the right-hand side of (3.23), we obtain 





d CF — G(é - 
oe cos} (€ + n) (é "|. (3.25) 
ds 2 
From (3.25), we find that except for a constant 
20 = F(é + n) — G(é — »). (3.26) 


By use of (3.14), (3.15), (3.16), and (3.26), the magnitude of the velocity w and 
the density p may be shown to be expressible solely in terms of w and the Mach 
number of the flow. Thus, from (3.2), (3.12) 


da OB da OB 
w= 1(= 45), v=/[_—+—}. (3.27) 
Ox Ox Oy doy 
Hence, by interchanging the independent and dependent variables, it follows that 
l fa oy lL fdx ax 
v= - (2-2), = —(=-=), (3.28) 
D\0B da D\0a_=saB 
where 
dx 0 dx 0 
i a es (3.29) 
da 0B OB da 


Computing the partial derivatives by use of (3.15), (3.16) and simplifying by use of 
(3.14), (3.26), we obtain 


u = 21 sin }w cos 6/sin w, v = 21 sin }w sin 6/sin w. (3.30) 
Hence, for the magnitude w of the velocity, we find 


21 sin }w 
———_: (3.31) 


w= 
sin w 


From (2,6), we see that 2 is equal to wj—a?. Making this substitution in (3.31) and 
dividing the resulting equation by ~, we obtain 

w 29/1 — (a;/w,)? sin fw 

eat / ( 1/ 1) 2 . (3.32) 


WwW sin w 





In Fig. 1, these curves are plotted for the following values of the Mach number, 
w;/a;=1.5, 2.0, 2.5, 3.0, 4.0. Note, by the discussion following equation (2.10), as 
w;/a, varies from 1 to ©, the permissible values of the upper bound of w/w; varies 
from 1.7 to 1. In each case, the upper bounds are to be determined by use of (2.10). 
The dotted lines in Fig. 1 denote these upper bounds. Other useful results may be 
obtained by combining (3.31) with (2.6) and (2.3). Thus, dividing the equation 
w*—a?=[ by w? and inserting the value of w as determined from (3.31) into the 
right-hand member of the resulting equation, we find after a few trigonometric sub- 
stitutions 





(3.33) 
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This equation is of value in determining the lower limit of the ratio w/a, namely, 
V2, for w=7/2. Again, inserting the value of w as determined from (3.31) into the 
. - s . . . . . 9 . 
left-hand side of (2.3), dividing the resulting equation by aj, and replacing the term 

o 2 2/2 > = 

P/at by wi/aij—1 (see 2.6), we obtain 
p sin w 
; . (3.34) 
pi 2 sin? $w./(w;/a)? — 1 

These curves are plotted in Fig. 2 for the values of the Mach number as indicated 


above. 
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An important case in practice’ (the jet problem) is that for which one stream line 
is a straight line. In this case, if we assume that the stream line coincides with the 
x-axis and is 7=0, (3.26) furnishes the result 


w 
wm 
Jt 


F(t) = G(é). (3.3: 


Under properly given Dirichlet data, the function F(&) can be determined and the 
representation of all stream lines can be obtained from (3.19), (3.20). Thus, if 0(&, c) 
is known along some known stream line »=c, then the equations (3.26) and (3.35) 


furnish the result 
20(é, c) = F(E+ c) — F(E — cc). (3.36) 


The equation (3.36) can be solved for F(£) by use of the theory of difference equations. 


’ J. Ackeret, Gasdynamik, Handbuch der Physik, vol. 7, pp. 318-322. 
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Unfortunately, 6(£, c) is unknown; 0(s), where s is the arc length parameter, is known. 
Hence, one must solve problems by an indirect method. That is, one must introduce 
a function F(£) and then determine the corresponding stream lines. 

In calculating the stream lines for particular functions F(&), the following analyti- 


cal-geometrical scheme appears to be the most satisfactory. First, obtain two curves 





7.0 








6.0 


























hk 
oO 





>| 




















2.0 


























of 





Fic. 2. 


of the generating Tschebyscheff net (one of each family) by use of equations (3.15), 
(3.16). That is, determine the curves 


xy -f cos F(t)dt, v1 -{ sin F(t)dt, (3.37) 


B 
Xe -f cos F(t)dt, Vo 


To obtain the initial point of each curve in the x, y plane, we compute one set of 
values of (x, y) by use of (3.15), (3.16), for some set of values of (a, 8) such as a=0, 
8=0. By translating the curves along each other, the complete Tschebyscheff net 
may be obtained. However, the translation must furnish curves of the families which 
correspond to equi-distant values of @ and B in order that the diagonal curves be 
stream lines. In view of (3.15), (3.16), this means that the abscissas of the initial 


points of two corresponding curves must be equal. 


3 
-f sin F(t)dt. (3.38) 
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Finally, as an example of this method, let us consider the case F(¢) =arc cos ¢. From 
(3.37), (3.38), it follows that 





a? i ee 1 
1=—) y= —“s — art Fare sin a, (3.39) 
6? B 1 
2_=—») = ——v/1 — BP? -—- — in p. 3.40 
v2 - V2 r V/1- 8 r arc sin 8 ( ) 





Fic. 3. 


The Tschebyscheff net and the resulting stream lines, obtained by the procedure 
outlined in the preceding paragraph, are illustrated in Fig. 3. By use of a protractor 
and the graphs of Figs. 1 and 2, the values of w/w; and p/p; can be immediately de- 
termined at each point of the plane. 

In concluding, it should be pointed out that it would be highly desirable to obtain 
a mechanical method for constructing the Tschebyscheff net when two stream lines 
are known. This would furnish a direct solution to the problem of the uniform flow 
of a supersonic jet. It appears that a more thorough understanding of the relation 
between Tschebyscheff nets and their diagonal curves is needed. 
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APPENDIX 


First, we shall give an analytic derivation of relations (3.6), (3.7). Since the pro- 
jections of the velocity vector on the normal to the characteristics have the magni- 
tude a, it follows that 


aA =>u——-—v—» (A.1) 
0a Oa 
fs] Ox 

wai ell wy eee sa ee (A.2) 
op 0p 


Further, by projecting the velocity vector on the tangents to the characteristics, we 


obtain 
Ox oy 











/wv2—v—aA=u—+r—, (A.3) 
Oa Oa 
Pe peg enn a (A.4) 
»: apap 
Solving relations (A.1), (A.2) and (A.3), (A.4) for u, v, we find 
0 0 J 2 S a. an 0 fa] 
w= - (a= 442) -% a B “(42 - 82), (A.5) 
D 0a op D op 0a 
ov oy 2 , - oF re] Ox 
p= — “(5 ae ~~) a (s=- A =), (A.6) 
D\ da aB D da ap 


where D is the Jacobian of the transformation (3.4). By interchanging dependent 
and independent variables, (A.5), (A.6) become 


w= 0(4— Bo) = VF e= a4 + p=), (A.7) 

Oy Oy Ox Ox 

v= ~ (4-8) - varpea aa + p=). (A.8) 
Ox Ox Oy dy 


From (A.7), (A.8) and (3.2), we find (3.6), (3.7). 

Next, we shall show that the Kérmén-Tsien relation (2.2) is valid, if the net of char- 
acteristics form a Tschebyscheff net. 

Since the net of characteristics is a Tschebyscheff net, we may consider the metric 
coefficients A and B as having the value unity. Hence, from (3.6), (3.7) and the chain 
rule for differentiation, we obtain 











oy O(a ~— B) oy O(a _ B) 

ee hen : —= ————-; A.9 
Ox e Ox oy ” oy , 
0 Scientia te) te) a ee 

2 fFEFAB O(a + B) | J/e+e—a %e +A) | (A. 10) 
Ox Ox oy oy 
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With the aid of (3.1), (3.2), we may write the integrability conditions of (A.9), 
(A.10) in the form 

fos Nees. 
—a/u? + v? — a? = 0, (A.12) 
On 


fe) 

— pa = Q, (A.11) 

OSs 
where 0/d0s represents differentiation along a stream line and 0/0n represents differ- 
ential along an equipotential. If we assume that a relation exists between p and p, 
then by use of the definition of a? (defined as dp/dp), we find from (11) 


i? d 
Pl ae hs (A. 13) 
dp’ dp 


Further, from the generalized Bernoulli relation 


ap _ 1 dp vm 


u?> + vy? + 2P(p) = constant, = 
dp p dp 
we find that (A.12) reduces to (A.13). Integrating (A.13), we obtain (2.2). 

The author wishes to thank Professor W. Prager and Dr. L. Bers of Brown Uni- 
versity for valuable criticisms and suggestions. Further, he is indebted to the Re- 
search Institute of the University of Texas for a grant which permitted the construc- 


tion of the diagrams. 
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ON THE STABILITY OF TWO-DIMENSIONAL PARALLEL FLOWS 
PART I.—GENERAL THEORY* 


BY 


C..<.. 45 
Guggenheim Laboratory, California Institute of Technology 


1. Introduction. The study of the stability of laminar motion and its transition 
to turbulence dates back to the time of Helmholtz and Reynolds [46], and had 
already attracted great attention at the end of the last century.** Since that time, 
the subject has not only become a major problem for workers in hydrodynamics, 
but has also attracted the attention of people like Lord Rayleigh [43-45], Lord 
Kelvin [20-21], Lorentz [29], Sommerfeld [58], and Heisenberg [14], whose chief 
interest is not limited to the study of mechanics. Although numerous contributions 
have since been made, the subject has remained one of considerable dispute, as can 
be seen from the two general lectures given by Taylor [70] and by Synge [63] as 
late as 1938. Still more recently, there appeared the work of Gértler [8, 9] and of 
Thomas [71]. 

Most of the work on the stability of laminar motions has the following final aims. 

1) The first aim is to determine whether a given flow (or a given class of flows) 
is ultimately unstable for sufficiently large Reynolds numbers. For this purpose, it 
is desirable to obtain some simple general criterion which will give a rapid classifica- 
tion of velocity profiles according to their stability. 

2) The second purpose is to determine the minimum critical Reynolds number at 
which instability begins. It is often easier to find sufficient conditions for stability 
than to find the condition for passage from stability to instability. 

3) Finally, we want to understand the physical mechanism underlying the 
phenomena by giving theoretical interpretations and experimental confirmations of 
the results obtained from mathematical analysis. 

Although numerous attempts have been made in these directions, especially for 
the apparently simplest cases of parallel flows in two dimensions, our knowledge is 
still very meagre. The classical case of plane Pouiseuilla motion has remained an un- 
settled problem,t and no satisfactory general results have been reached regarding the 
stability of a real fluid. The best-known general criterion is that of Rayleigh (1880) 
and Tollmien [74], classifying profiles according to the occurrence of a flextt with re- 
spect to the stability of a fluid at infinite Reynolds numbers. However, the sig- 
nificance of their results has been too much exaggerated and often misunderstood, 
and no physical interpretation has ever been given. The present work offers such an 
interpretation, but also shows that the results can only give some indication regarding 


* Received March 3, 1945. An abstract of this paper has already appeared under the same title [27]. 

** In 1888, the problem was proposed by Rayleigh and Stokes as the subject for the Adams Prize 
Essay. Cf. p. 321 of Ref. [21], and also the footnote on p. 267 of Ref. [44]. 

t Cf. Synge’s lecture [63]. 

tt Following Professors Frank Morley and H. Bateman, we shall use the word “flex” for “point of 


inflection.” 
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the instability of a real (viscous) fluid. This will be discussed in more detail below. 

The chief aim of the present work is to try to answer the three questions men- 
tioned above for two-dimensional parallel flows. This work is divided into three parts. 
Part I (the present paper) deals with the general mathematical theory, with par- 
ticular emphasis on attempting to clarify the mathematical difficulties involved in 
the solution of the equation of stability. Part II deals with the stability problem in an 
inviscid fluid (infinite Reynolds numbers). Part III deals with the problem in a real 
fluid. The following results have been obtained. 

1) It is shown that all velocity distributions of the symmetrical type and of the 
boundary-layer type are unstable for sufficiently large (but finite) values of the 
Reynolds number (Part III). The plane Poiseuille motion is included as a special 
case. 
2) A simple approximate method is obtained by which one can calculate the 
minimum Reynolds number marking the beginning of instability with very little 
numerical labor (Part ITI). 

3) The tendency toward instability of a profile with a point of inflection is inter- 
preted by considering the distribution of vorticity (Part II). The effect of viscosity is 
considered as diffusing the disturbance from the “critical layer” inside the fluid and 
from the solid boundary. A very simple quantity is thereby derived which serves as 
a measure of the effect of viscosity (Part III). This can also be easily connected with 
the general mathematical theory. 

As numerical examples, we have worked out the curve of neutral stability for the 
Poiseuille case and the Blasius case. Comparisons with existing results are discussed 
(Part III). The relation between instability and transition to turbulence is also dis- 
cussed in Part III of this work. 

Since some of the present results differ markedly from customary beliefs, it is 
necessary to trace the history of the existing lines of thought in order to give proper 
recognition to earlier ideas and results used in the present work, and to analyze all 
the results in disagreement with present conclusions. This requires the repetition of 
some known results when they fall into the present line of treatment. The review of 
literature is not intended to be exhaustive; only the necessary references are cited. 
A more complete bibliography up to 1932 has been given by Bateman [2]. 

2. Historical survey of existing theories. There seem to be two schools of thought 
in regard to the cause of transition from steady to turbulent conditions. One school 
contends that transition is due to a definite instability of the flow, i.e., to a condition 
in which infinitesimal disturbances grow exponentially. The second school regards 
the motion in most cases as definitely stable for infinitesimal disturbances but liable 
to be made turbulent by suitable disturbances of finite magnitude or by a large 
enough pressure gradient. Both schools, however, generally agree that the fluid can 
be considered as incompressible and that its motion is governed by the Navier- 
Stokes equations of motion. Since the agreement between theory and experiment has 
not been very satisfactory, it has also been proposed that the cause of transition 
must be traced back to the effect of compressibility or to the possible failure of the 
Navier-Stokes equations. The present work tends to confirm the simplest point of 
view that the motion in most cases is definitely unstable for infinitesimal disturbances 
governed by the Navier-Stokes equations for an incompressible fluid. 

The theory of finite disturbances dates back to Reynolds [46] and Kelvin [21]. 
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It was developed by Schiller, Taylor and others.* Mathematical investigations of 
such finite disturbances are mainly based on considerations of energy or of the square 
of the vorticity of the disturbance, because the solution of the non-linear equations 
satisfied by the disturbance is extremely difficult. At the end of Part III we shall 
briefly discuss this line of thought together with the results of the present paper. For 
more details, the reader is referred to the lecture of Taylor [70] and the papers of 
Synge [62] and Thomas [71]. 

For small disturbances, positive definite integrals of the energy and vorticity of 
the disturbance have been extensively used. These considerations have been discussed 
by Orr [37], Lorentz [29], von K4rm4n [18], Synge [63, 64] and others. For excel- 
lent accounts of this phase of the theory, the reader is referred to the works of 
Noether [35], von Kérm4n [18], Prandtl [42], and Synge [64]. Additional references 
are cited at the end of this paper. As is now well-known, this method can only give 
sufficient conditions for stability. Also, since all disturbances are usually allowed, in- 
cluding those which do not satisfy the hydrodynamic equations of motion, a larger 
viscous decay is required to insure stability than when these disturbances are ex- 
cluded. Consequently, the limit of stability is always found to be much lower than 
that indicated by experiment. However, from these considerations, Synge [63] has 
arrived at a very convenient form of a sufficient condition for the stability of two- 
dimensional parallel flows with respect to two-dimensional disturbances. This will be 
found very useful for the discussions in Part III. 

To get more concrete results, we have to solve the linearized equations satisfied 
by the disturbance. The most successful case appeared to be Taylor’s treatment of 
Couette flow [67] between concentric cylinders. His work was verified by the experi- 
ments carried out by himself [67, 69] and by others [28]. A rigorous mathematical 
investigation in this connection was made by Faxon [4]. In fact, it is now known 
that his analysis is a typical case of the stability of a fluid motion where the centri- 
fugal force plays a dominant part. Such cases were first considered by Lord Rayleigh 
[45], who gave a condition for the stability of an inviscid fluid. Mathematical proof 
of a sufficient condition of stability of Couette flow was recently given by Synge [65]. 
Extension of Taylor’s work to the boundary layer over a curved wall was carried out 
by Gértler [8, 9], who used numerical methods successfully. 

While the investigation of curved flows was uneventful, the investigation of axi- 
ally symmetrical flows was not extensive. The Poiseuille flow in a circular pipe was 
studied by Sex! [55] with a conclusion of stability. Prandtl [42] gave some discus- 
sions of the possible cause of instability in his article in the book “Aerodynamic 
Theory,” edited by Durand. 

The most extensive discussion of hydrodynamic stability seems to be the treat- 
ment of parallel flows by attempting to solve the eigen-value problem associated 
with the linearized equations governing the disturbance. This line of development 
can be easily traced in the work of Helmholtz, Lord Rayleigh [43, 44], Orr [37], 
Sommerfeld [58], von Mises [31, 32], Hopf [16], Prandtl [41], Tietjens [72], 
Heisenberg [14], Tollmien [73-75], and Schlichting [52-54]. Other contributions are 
those of Noether [36], Solberg [57], Southwell [59], Squire [60], Goldstein [6], 
Pekeris [39, 40], Synge [61-65] and Langer [25]. 


* See Taylor's lecture [70] for references to the works on finite disturbances. 
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For convenience, the theory deals with two-dimensional wavy disturbances propa- 
gated along the direction of the main flow. Squire [60] has shown that three-dimen- 
sional wavy disturbances are more stable than two-dimensional ones. However, 
Prandtl still mentions the possibility of greater instability of three-dimensional dis- 
turbances in his article [42] appearing after Squire’s paper. 

The first study of two-dimensional hydrodynamic stability seems to have been 
made by Helmholtz. He proved the instability of wavy disturbances over the surface 
of discontinuity of two parallel streams of different velocities. Later, Rayleigh [43 | 
realized that Helmholtz’s approximation was not good enough to bring out the main 
features of a flow with continuous velocity distributions. He therefore made an im- 

proved approximation consisting 

















(a) (b) (c) of several linear profiles joined up 
+t i continuously. The vorticity dis- 
. tribution then has constant values 
pI in several layers, but has a dis- 
continuity in passing from one 
; : ; layer to another. Investigations 
, | t 4 with continuous vorticity distri- 
' | a ° ‘ ; 
1 oi ny butions were also made. Rayleigh’s 








work was mainly concerned with 
Fic. 1. Broken profiles investigated by Lord Rayleigh. Case an inviscid fluid. Two main results 
(a) may be unstable; the other two cases are stable. were obtained. The first is that in- 
stability (in an inviscid fluid) can 

only occur with velocity distributions having a point of inflection. It is usually be- 
lieved that Lord Rayleigh has also proved that damped disturbances can also only 
occur with such profiles. The possibility of a disturbance for a profile without a flex 
then becomes a paradox [5]. Actually, Rayleigh’s proof does not lead to such a 
conclusion. This point will be more fully discussed in §5 and Part II. Rayleigh’s 
second result is obtained from the analysis of broken linear profiles; it substantiates 
the first result by demonstrating definite instability of broken linear velocity distri- 
butions of the type shown in Fig. 1(a), and only stability in the other cases. Rayleigh 
[43] supported his result by obtaining the condition determining stability in the 


f (w — c)-*dy = 0, (2.1) 


yl 


approximate form 


where w(y) is the velocity distribution, y; and ye are the coordinates of the solid 
boundaries, and c is a constant the real part of which represents the wave velocity 
and the imaginary part of which gives damping or amplification. 

Meanwhile, the exact analysis of linear velocity distributions including the effect 
of viscosity was given by von Mises [31, 32], and Hopf [16] and was also studied by 
Rayleigh [44]. The results indicate only stability. Prandtl and Tietjens [72] applied 
Rayleigh’s method of approximation to the stability of the boundary layer, taking 
account of the effect of viscosity. In such an approximation, the inner friction layer 
mentioned above ($1) for continuous vorticity distributions is left out. The result of 
Tietjens did not give a minimum critical Reynolds number. 

It was Heisenberg [14] who first successfully studied the stability of a variable 
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continuous vorticity distribution. As a particular example, he demonstrated that the 
plane Poiseuille flow was unstable for sufficiently large Reynolds numbers. Also, using 
the same equation (2.1) with which Rayleigh supported his approximation with linear 
profiles, Heisenberg pointed out the fallacy in Rayleigh’s method. The essential point 
is that the corners in the velocity profile introduce extraneous roots of the above 
equation for c. Consequently, the results of this type of analysis depend upon the 
manner in which the velocity distribution is approximated. 

Heisenberg’s numerical computation was, however, incomplete and very rough, 
and his theory was not generally accepted. Better known are the results of Tollmien 
and Schlichting. They studied the cases of Blasius [73] and plane Couette flow [48], 
using Heisenberg’s theory essentially. The former case was pursued very much in 
detail. For the latter case, Schlichting followed the idea of Prandtl, asserting that the 
instability may be attributed to the initial unsteady distribution prior to the forma- 
tion of the linear profile. Indeed, the same kind of idea was also suggested by Prandtl 
to account for the instability of Poiseuille flow by ascribing it to the entrance section 
where the profile is not yet parabolic [41]. This problem will be discussed in some de- 
tail later (§14, Part ITI). 

For an inviscid fluid, Tollmien has also proved the instability of boundary-layer 
and symmetrical profiles with a point of inflection [74]. For a viscous fluid, the pres- 
ent investigation shows that instability depends upon the general type of these pro- 
files rather than on the appearance of the point of inflection. The inner friction layer 
plays a dominant role in determining the instability. Attempts to interpret this point 
physically are given by Prandtl [42] and in the present paper. 

3. General formulation of the problem. We shall now formulate the problem of 
the stability of two-dimensional parallel flows mathematically. In the first place, we 
note that if the steady motion is strictly two-dimensional and parallel, the velocity 
distribution must be either linear or parabolic (if body forces are absent). We then 
have one of the following: 1) the plane Couette flows; 2) the plane Poiseuille flow; 
3) a combination of these two flows. The problem would then be very restricted. 

However, there are a large number of cases where the flow is essentially parallel 
to one direction. These are the cases where the boundary-layer consideration is per- 
missible. The following are important special cases belonging to this class: 4) inlet 
flow between parallel walls, flow in a slightly convergent or divergent channel; 5) flow 
along a flat plate; 6) wake behind a cylindrical body, jet from a narrow slit. Whether 
these flows can be properly considered as belonging to the same class as the above 
three is a question of some controversy. Taylor has criticized Tollmien’s work with 
the boundary layer on this ground [70]. In the Appendix to Part III of this work, 
we shall try to demonstrate that this treatment is generally permissible, but that the 
interpretation of the results must be taken up with care. A discussion of Tollmien’s 
work will also be found there. 

In considering the stability of the main flow, we superpose upon it a hydrody- 
namically possible small disturbance, and consider its behavior. The disturbance is 
small in the sense that the inertia forces corresponding to the disturbance alone are 
negligible and that its behavior is unaltered when its amplitude is (say) doubled or 
halved. It is then simplest to consider separate harmonic components with respect to 
time, which may be damped, neutral, or self-excited. By considering disturbances 
which are also spacially periodic both in the direction of flow and in the direction 
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perpendicular to the plane of symmetry of the main motion, Squire [60] was able to 
show that two-dimensional disturbances are less stable than three-dimensional dis- 
turbances. Hence, important features of the stability problem can be obtained by 
considering two-dimensional disturbances alone. This is an essential difference be- 
tween the stability of a parallel flow and of a curved flow. In the latter case, three- 
dimensional disturbances are of utmost importance. 

The consideration of periodic disturbances alone is again a question of some con- 
troversy. Justification has been attempted and objection has been raised. We shall 
see later that at least the difficulties raised are chiefly caused by a misinterpretation 
of the mathematical results. 

Admitting that we can consider two-dimensional disturbances alone, we have a 
much simplified physical picture at hand. If the effect of viscosity is negligible, we 
have the well-known fact of conservation of vorticity for two-dimensional motions. 
Actually, the stability problem is found to depend both on the inertia forces and on 
the viscous forces. However, the effect of viscosity is also well-known to be one of 
diffusion of vorticity. Thus, important results can be expected from considerations 
of vorticity transfer. 

Let us now proceed with the mathematical formulation of the problem. We shall 
give a complete derivation of the stability equations so that we can see how to settle 
the disputes about the approximations in considering velocity distributions of the 
boundary-layer type. 

Admitting Squire’s work as a proper indication that only two-dimensional dis- 
turbances need be considered, we may conveniently consider the equation of vorticity 


Avs + VuAve — ved, = vAAY, (3.1) 
with the velocity components 
Oy Oy 
u=Wy=—) p= —ye= =, (3.2) 
o% Ox 
and the vorticity 
ov Ou 
eS =O = Cee T Pe ~~ BP. (3.3) 
Ox oy 


As usual, v is the kinematical viscosity. We may add that Squire’s original proof was 
intended for flow bounded between two parallel walls. There is no difficulty in seeing 
that the proof holds also for a fluid extending to infinity,* because the boundary con- 
ditions for the disturbance are essentially the same. 

Let us put 

y = V(x, y) + ¥’(x, y, 4, (3.4) 
where V(x, y) represents the steady main flow and w’(x, y, ¢) represents the dis- 
turbance. Main flows which vary but slowly with time can also be treated this way, 
but we shall restrict ourselves to steady flows in order to fix our ideas. 

If we substitute (3.4) into (3.1) the terms corresponding to the main flow cancel 
out. If we then drop the terms quadratic in y’(x, y, 4) and its derivatives, we have 
the equation 

Ayi + W,Ayz — VApy + WAV, — YZ AV, = vAdy’. (3.5) 


* Cf. Ref. [15] 
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We shall now assume the flow to be essentially parallel to the x-axis. Using the bound- 
ary-layer approximation, we should drop the x-derivative of any quantity connected 
with the main flow compared with its y-derivative. But for the disturbance we would 
expect ¥/ and W/ to be of the same order of magnitude. This will be verified a 
posteriori in the specific examples. Further discussions will be found in the Appendix 
to Part III. With these considerations, (3.5) reduces to 


3 
Ayi + V,Ay/ —yvZ — vdAy’. (3.6) 
V 
Now we shall make an approximation of the same order by taking for w=WV, and 
0°w/dy? =0°V /dy® their local values at a given value x» of x. Then we may write 


Ayi + w(y)dys — w'(y)ps = vAdy’. (3.7) 


For the boundary conditions, we shall also consider the /ocal boundaries. The problem 
is then essentially simplified, and can be treated similarly to plane Couette and 
Poiseuille flows. We consider a main flow between two parallel planes y= y, and y = y2 
with a more or less arbitrary distribution of velocity w(y). Then the disturbance 
W(x, y, t) must be found as a solution of (3.7) satisfying the conditions u’ =v’ =0 
over the boundaries. 

The usual way of dealing with the solution of (3.7) subject to given boundary con- 
ditions is to consider periodic disturbances. We shall refer all velocities to a charac- 
teristic velocity U and all lengths to a characteristic length /, and define the Reynolds 
number R = Ul/v. The two-dimensional periodic disturbance of a field of flow in which 
the main flow is w(y) may be represented by the stream function ~’ =¢(y)e‘#(=-«, 
and the linearized differential equation for ¢(y) is 


 . 
(w — c)(¢” — a®d) — wb = — — (¢'* — 20°” + a9), e (3.8) 
aR 
as can be easily obtained from (3.7). We shall take a always real and positive, while c 
may be complex; thus, 


C= 6 + 1G. (3.9) 


To fix our ideas about the boundary conditions, let us consider a flow between the 
planes y=y; and y=ye. The equation (3.8) is then to be solved under the boundary 


conditions 
(v1) = 90, o(y2) = 9, O(n) = 9, = G'(y2) = 0. (3.10) 


Let us now forget about the physical problem and consider the differential equa- 
tion (3.8) as a linear differential equation of the fourth order in the complex y-plane. 
To be sure, the function w(y) is defined only for real values of y between y; and ye. 
We can of course, consider it as defined for other values of y by analytical continua- 
tion. We shall assume that the function thus defined is holomorphic in every finite 
region with which we shall be concerned. The equation (3.8) then has every point in 
the region under consideration as a regular point, and its coefficients are also entire 
functions of the parameters c, a, and aR (regarded as complex variables). By a well- 
known theorem in the theory of differential equations, there exists a fundamental 
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system of four solutions of (3.8) which are analytic functions of the variable y and 
of the parameters c, a, and aR, being in fact entire functions of the parameters. The 
consequences of these simple general analytical considerations appear to have escaped 
serious attention from earlier investigators. In §§4, 5 of this paper, we shall find this 
type of consideration very important in settling the controversies about the question 
of convergence of the series used in the actual solution of equations (3.8) and (3.14). 

Let us denote the above-mentioned system of solutions of (3.8) by ¢:(y), d2(y), 
o3(y), d4(y), the dependence upon the parameters c, a, aR being understood. The con- 
ditions (3.10) then give rise to the secular equation 


O1(¥1) ely) slvr) = hays) 

o1(¥2) $2( V2) o3( V2) o4( 2) : ' 
. = 0. (3.11) 

oi (¥1) 2 (¥1) O38 (91) O64 (1) 

di (ye) o2(¥2) 3 (v2) 4 (2) 


—_ 
— 


Fic, a, aR) = 


Since the function F(c, a, aR) is an entire function of the variables c, a, aR, we may 


solve for c, obtaining : 

c = cla, R). (3.12) 
There may be several branches of the solution, or there may be none as in the case 
when F(c, a, aR) is (say) exp(aRc). In general, we would expect the solution to be 
unique, or we may consider only one branch of the solution. 


Since a and 2 are later taken to be real and positive, it is convenient to separate 


(3.12) into its real and imaginary parts. Thus, 


Cr = ¢-(a, R), c; = ci(a, R). (3.13) 


It is customary to plot curves of constant c; or ac; in the aR-plane. The curve c;=0 
gives the limit of stability. 

Weare particularly interested in the case when the Reynolds number is very large. 
The study of this case is complicated by the fact that the functions $1, deo, $3, $4 in- 
volved have essential singularities at the infinite point of the aR-plane. From the 
differential equation (3.8) itself, we see that when aR— ~, we have the equation 


(w — c)(o” — a*d) — wo = 0, (3.14) 
which is only of the second order. Thus, two solutions of (3.8) are lost. From detailed 
mathematical investigations, we shall find later that two linearly independent solu- 
tions of (3.8), say ¢: and ¢2, will satisfy (3.14) in the limit of infinite aR, except 
along certain straight lines through the point w=c. The other two linearly independ- 
ent solutions ¢3 and ¢, are highly oscillating for large aR and would therefore disap- 
apear in the limit of infinite aR. Furthermore, we shall see that ¢3 and ¢, can be so 
chosen that if }3(y1)>@4(y1), then $3(v2)<.(y2), with corresponding relations for their 
derivatives. It then appears plausible that the limiting form of (3.11) for infinite aR is 


| o1(¥2) $2(91) 
| p1(y2) $2(y2) 


(3.15) 


j= 
| 


with ¢:(y), 2(y) satisfying (3.14). 
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The condition (3.15) states that we are looking for a solution ¢(y) of (3.14) satisfy- 
o(y1) = 0, o(y2) = 0, (3.16) 


with the other two conditions of (3.10) relaxed. Physically, this means that we allow 
a slipping along the walls y=y, and y=ye2. For very large Reynolds numbers, only a 
very thin layer of fluid will stick to the solid, and we have naturally an apparent 
slipping. These points will be taken up again more carefully (§6) after a thorough 
mathematical investigation of the solutions. 

4. Solution of the equation of Orr and Sommerfeld by methods of successive ap- 
proximation. The stability equation of Orr and Sommerfeld 


i 
(w — c)(o” — a*p) — wd = — a: (p'’ — 2a*p”’ + ap) (4.1) 
a 


has a fundamental system of four solutions, which are analytic functions of y (where- 
ever w(y) is analytic) and which are entire functions of a, c, and aR. In order to obtain 
useful solutions, it is usual to expand the solutions as power series of a suitable small 
parameter, say, (aR)~!. However, since (aR)~! occurs with the highest derivative 
in (4.1), the study of such an expansion becomes very complicated. It will be done 
later. 

a) Solution by convergent series. An alternative method* is to choose a small pa- 
rameter € related to (a2R)~' and first make a change of variable (yo being an arbi- 


trary point so far) 
y— You en, o(y) = x(n), (4.2) 


so that (4.1) becomes 





(w — c)(x” — are*x) — ew"’x = — (xi¥ — 2ae2x”’ + atetx), (4.3) 
aRe* 


where 
, 0 2 
w— ¢ = (Wo — c) + wo (en) Sa ee ts , 
pan (4.4) 


eS 
I 


WwW 
” = we’ + wo’ (en) ae OP tt 0 


The solution is then obtained in the form 
o(y) = x(n) = x(n) + exM(n) + OxO(n) +---, (4.5) 


and the differential equations for the approximations of successive orders can be ob- 
tained by substituting (4.4) and (4.5) into (4.3) and equating all the coefficients of 


the various powers of € to zero. 
If we take yo to be the point where w=c, the proper choice of the parameter € is 


¢ = (aR)~!/3, (4.6) 

The differential equations for the functions x(n), x(n), x(n), - - - areas follows: 
0. we (0/7 + iy (iv = 0, 

wlan" + | és 


ewe x” + ixMiv = La a(x), (un 2 1), J 


* This method was first used by Heisenberg, loc. cit. [14], p. 588. 
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where L,~1(x) is a linear combination of x(n), x(n), - - - , x°"7(n) and their de- 
rivatives. In particular, 

Lo(x) = wo’ (x — 32x"), - (4.8) 

We note that the homogeneous part is the same for all the differential equations 

in (4.7). Hence, if we can solve for the first approximation, the rest can all be obtained 

by guadratures. Indeed, the first equation of (4.7) is Stokes’ equation* for x’, and 

its solution can be readily expressed in terms of Bessel functions of the order 1/3. 


Thus, for the first equation of (4.7) we have the four particular integrals** 


\ 


n ” 
0 (0) 1/2 Gita ss 3/2 | 
X1 = Mn, X3 -{ an f dnn Hy, 3| 3 (iaon) |, | 
Loo +00 | 


(4.9) 
0 (0) ” " 2 c054 , 3/2 
x2 =1, “i = dn dnn Hy 3 [3 (iaon) |, 
—o —o 


/ /3 
ao = (Wo ye “ 


where 
(4.10) 


The higher approximations are given by 


(n) 7 7 7 (O)rs ’ Ors Oss ’ (O)y7 
x = dn dn 4x4 dnxs  Ln-1(x) — xs dnxs Ln-i(x)?, 
6 (4.11) 


(é = 1,2, 3,6). 


These are the explicit formulae for finding the approximations of various orders. In 
actual calculations, only the initial approximation (4.9) is required. Furthermore, the 
series (4.5) is convergent provided ¢ is restricted so that the series (4.4) are convergent. 
For then the differential equation (4.3) for x(n), when normalized, has analytic func- 
tions of the parameter € as its coefficients. Hence, a fundamental system of its solu- 
tions consists of four analytic functions of e. 

It should be mentioned that if yo is not taken at the particular point for which 
w=c, the proper parameter to be taken is (aR)? instead of (aR)—"*. In this case, 
all the approximations can be expressed in terms of elementary transcendental func- 
tions. However, it is not found particularly advantageous to do so, because the study 
of “crossing substitution” (§5) would not be easy. Also, the method is then too much 
different from those used by earlier investigators to allow an easy comparison of 


the results. 
b) Solution by asymptotic series. Although the previous method is theoretically 





* Cf. the exact treatment of (4.1) by Hopf [16] and Rayleigh [44] for the case w’’ =0. 

** Note that x,” and Xs and also x; have no branch point at 7=0. The order of the solutions 
| b1, b2, os, 64} agrees with Tollmien’s notation. They are | 3, ¢4, ¢1, ¢2} in Heisenberg’s notation. Heisen- 
berg gave the solutions ¢; and ¢,4 in terms of Hankel functions in the form 


1 2 
$j; = (w — Cc) f- ws aa) fn, (j = 1, 2), 
9 3 


(p. 289, and Eq. (19a) p. 591). It can be easily verified that these are the same as ne up to a constant 
factor wo’e and to the proper order of approximation. Note that throughout Heisenberg’s paper, 7 is to 
be replaced by —iin order to conform to our notation. This can be seen from a comparison of our Eq. (4.1) 
with his Eq. (7a). The difference arises from a difference of notation in the stream function y'(x, y, ¢). 
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complete, it is usually more convenient to use asymptotic series for numerical pur- 

poses, particularly in dealing with boundary-value problems. Heisenberg has given 

two asymptotic methods, each of which gives only two particular solutions of (4.1). 

These methods will now be described and investigated mathematically in more de- 

tail, because Heisenberg’s work has received some criticism in this connection.* 
The first of these methods is to develop ¢(y) in powers of (aR)~'. We put 


o(y) = G6 (y) + (aR) (y) + (aR)-*'(y) tes, (4.12) 
and substitute in (4.1). Comparing corresponding powers of (a#R)~', we have the fol- 
lowing differential equations 
(w— c)(p” — a®g™) — wo = 0, 

(w — c)(o®” — ad) — wo = — i [po e-Yiv — 2a2p kD" + aig (td), (4.13) 
(k 2 1). ) 

The initial approximation satisfies the inviscid equation and can be solved by de- 

veloping ¢ in powers of a®. Indeed, two particular integrals of (4.13) are 


(0) 


(w — c)[ho(y) + a hoy) + a hs(y) +--+], 





79) = 
z : P (4.14) 
» = (w—olkily) ta ks(y) +a ks(y) +--- J, J 
where 
ho(y) = 1, htons2(¥) = i dy(w — of dy(w — c)*hen(y), 
yl vl 
(n = 0), 
> (4.15) 
y 7] y 
ki(y) = dy(w — c)~*, Ren+s(y) -f dy(w — of dy(w — €)~*Rongi(y), 
“9 yl yl 
(n = 0). | 


The point y; might have been any fixed point instead of one of the end points; but 
it is found convenient to take it this way. 

Having found two particular integrals for ¢, we can obtain the higher approxi- 
mations by guadratures. In actual calculations, this is not necessary 

Because of the general nature of the Eq. (4.1), @(y) is an entire function of aR. 
Hence, the infinite point of the aR-plane is a singular point, unless ¢(y) is independ- 
ent of aR. Consequently, the series (4.12) is asymptotic, unless ¢(y) is a polynomial 
in (aR)~!. We note also that (4.13) is of the second order, so that only two solutions 
are obtained by this method. The solutions of (4.13) are entire functions of a? and 
hence the series (4.14) are uniformly convergent for any finite region of the complex 
a?-plane, for a fixed value of y, except when y is the singular point yo of the differential 
equation (4.13).f 

* Tollmien, loc. cit., 1929, p. 43. 

t This can also be seen from the series itself. So long as it is possible to run a path of finite length 


from y; to y on which w—c 0, the general terms a?"hz, and a***1k», 4; of the two series are bounded by 
A(aM)2"/(2n)! and B(aM)2"*!/(2n+1)!, respectively, (A, B, M being (suitably) fixed constants), and 
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In fact, the differential equation (4.13) has a logarithmic singularity at the 
point yo. This point is, however, an ordinary point in the exact equation (4.1), 
and the singularity is introduced purely by the method of asymptotic integration. 
However, the appearance of this singularity gives a serious ambiguity in the deter- 
mination of the correct path leading from y; to y in order that (4.14) may give valid 
approximations to integrals of (4.1) all along the path.* The proper way to settle this 
question is to compare the solutions (4.14) with the asymptotic expansions of the regu- 
lar solutions obtained by the previous method. This will be done later after we have 
described the second asymptotic method of Heisenberg for the other two particular 
integrals; for the same kind of problem also arises there. 

To obtain two other integrals of (4.1) in asymptotic forms, let us make the trans- 


formation 
o@ = exp { f sash. (4.16) 


Then, we obtain the non-linear differential equation 


me 


~ ‘ , 9) 
(w—c) (ge? + 2') —a*?} —w 


1 
= — — {gett 62g’ + 3g’2 + 4ge” + go’ — 2a*(g* + 2’) + a} (4.17) 


for the function Py We try to solve this by putting 
g(y) = (aR)!/*go(y) + girly) + (aR)~"/*g0(y) + °°. (4.18) 
Then, we obtain the set of equations 
2 4 aie 2 
(w — c)go = — igo, w— c)(go + 2gog1) = — t(4gogi + Ogoge ), 
w — ogi + gi + 2goge — a) — w = — i(dgoge + 6gi + 2gogigd + 3g0 — 2a g0), 


Hence we can obtain the successive approximations without integration. Thus, 


_— 5 go 
go=tvVil(w-—c), feos =n (4.19) 
2 £0 


For definiteness, we define 


hence the series converge like the cosine and the sine series, respectively. Heisenberg did not prove the 
convergence of these series, but stated that their convergence can be hoped to be sufficiently rapid for a? 
of the order of unity (loc. cit., 1924, pp. 584, 587). This was made a point of criticism by Tollmien (loc. 
cit., 1929, p. 43). 

* Considerable dispute has arisen in this connection. Note that it is impossible to dispense with this 
difficulty by remarking that the two different determinations will differ only by a constant multiple 
of a particular integral. If we draw two paths from y, to y and obtain such a difference in the solution, it 
is evident that the asymptotic solution cannot be valid on both paths, because the exact equation (4.1) 
has no singular point at y=yo and hence its solution must be single-valued. Although a mistake here 
would not cause serious difficulties so far as the numerical evaluation of the eigen-value problem is con- 
cerned, it does lead to misunderstanding and confusion elsewhere. Even after Heisenberg and Tollmien 
have analyzed this problem in some detail, they still take the very misleading step of taking the complex 
conjugate of the inviscid equation (Heisenberg, loc. cit., 1924, p. 596; Tollmien, loc. cit., 1935, p. 88). 
This point will be discussed more fully later. 
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T 
arg i = cs arg(w—c)>O for w—c>0. (4.20) 
For negative values of w—c, we cannot decide, without further investigation, whether 
arg (w—c)=+7 or arg (w—c)=—7. The point yo, where w=c, appeared in the 


previous asymptotic solution as a logarithmic branch point; here it-is an algebraic 
branch point. The determination of the correct path should follow the same criterion 
as the other two integrals, that (4.18) gives two asymptotic solutions of the exact 
equation (4.1) all along the path. This path might be expected to be the same as that 
in the previous case. All these will be discussed in the next section. 

After such a question is settled, substitution of (4.19) into (4.16) and (4.17) gives 


the two asymptotic solutions 


o3(y) = (w — c)~*/4 exp \_ f \/iaR(w — c) ay} ; 
. y 


| : } (4.21) 
oa(y) = (w — c)-*/4 exp +f ViaR(w — c) ay\ , | 
vo ) 


where factors of the order exp (aR)~”?=1+0{(aR)-"?} are taken as unity. 

5. Analytical properties of the solutions. Having thus obtained four asymptotic 
solutions of the equation (4.1), we must try to correlate them with the four solutions 
(4.9) and (4.11), and above all to study the correct determination of path around 
the artificial singularity introduced by the asymptotic methods. For this purpose, we 
consider the asymptotic expansions of the four regular solutions obtained by the first 
method and transfer back to the independent variable y. 

Let us recall that the asymptotic expansions of the Hankel functions HE(é), 
H)() are given by [76], 


_ 


a J 2 1/2 FP oT 
Hy3);3(4) ~ | — exp <21€ — — 
rt \ 12 





Lia FOO) 
b fh + \. 


r=1 (2i&)" 


| 
(—a <argéi < 2m), | 
a A : } (6.1) 
nto (2)"aw {4 2h 4 E824, | 
3(¢{) ~ [ - exp 4(—i{[ — — ——————->, | 
. oS ie 12/) a (288) | 
(—2n <argit<n). } 
If we put 3& = 2(taon)*/?, then (5.1) becomes 
(1) 3/2 3 a/2 3 3/2 5a 52 —3/2 
Hi [3 fags) —_ (=) (tan) “exp {i(a0n" e oo —\ }1 + O(n y}, 
TT ~-/ 
i” 53 
( = = < arg (aw oe ed by 
6 6 
‘ore 
5\ —3/2 
)}, 


(2) 3/2 sy =a 3/2 4i Sr 
A, 3[2(éexon)” }~ (=) (icon) : “exp {3(aen) sac + “ {1+ O(n 
T 


llr T 
(- — < arg (aan) < * . 
6 6 
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With the help of these formulae and using the legitimate process of integrating 
the asymptotic expansions term by term, we obtain 


ane! ant? \ 
(0) . (1) €Wo 2 1 P Wo . 
mn re Yet Te ——-1S teh 
2wo We € 2 
aati om 
(0 (1) wp “0 
xe texe ~1l+e—yq logy ~ 1 + — » log y, 
Wo Wo 
(0) —5/4 (2 3/2 Sat 4) 
x3 ™ const. 7 €xp 1 3(a0on) e f 
(5.3) 
5/4 y ee ee i] 
= const. (y — yo) exp {- f V iaRw¢ (y — yo) dye , | 
yl J 
(0) —5/4 iz , 3/2 wi/4) | 
X4 ™ const. 7 exp 13(aon) e { 
—5/4 ¥ = = a ————————— 
= const. (vy — Yo) exp f V iaRwy (vy — vo) dy? . 


vo 


These formulae can be easily seen to agree with the four asymptotic solutions (4.13) 
and (4.21) to the proper order of approximation, if we replace y; by yo in ¢{? (which 
is permissible). 

In evaluating the asymptotic expressions (5.3), the argument aon must satisfy 
both requirements specified in (5.2), i.e., 


— 7x/6 < arg (aon) < 7/6. (5.4) 


In this range, the asymptotic solutions (4.13) and (4.21) hold. Having thus established 
the range of validity of these solutions, we no longer need to make further com- 
parisons of the two methods of solution. 

At least three plans are now possible for further numerical work. First, we may 
use the four solutions obtained in the approximate form (4.9). Secondly, we may 
use the four asymptotic solutions (4.14) and (4.21). Thirdly, we may approximate 
ib1, be, 3, bs} by the four functions | ¢{”, o!, x}, xf} given by (4.14) and (4.9). 
The first method 1s very similar to the method used by Hopf [16] and Tietjens [72] 
for linear velocity distributions, where the exact solutions are given by functions of 
the general nature of those in (4.9). For curved velocity distributions, the functions 
x\”, xo” do not give ¢; and ¢» with sufficient accuracy, and this plan is not good. The 
second plan was used by Heisenberg in his investigation of the stability of the 
Poiseuile flow; but he also realized that it served only part of his purpose, and he 
stated that the third plan should be used.* Tollmien substantially adopted the third 
plan for his investigation of the stability of the boundary layer, although he did not 
point out the connection of his method with Heisenberg’s work. Instead of the expres- 
sions (4.14) for ¢: and ¢2, he used solutions in the forms of power series in y. These 
solutions are easily manageable only for linear and parabolic velocity distributions. 
Accordingly, he tried to approximate the Blasius profile with such profiles. Since such 
approximations are not good enough in the neighborhood of the point y= yo, where 
w=c, his discussion becomes very complicated. In the present work, we base our cal- 


* Loc. cit., p. 404. 
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culations upon the use of (4.14). It will be seen that our method can be applied to any 
profile with good accuracy. A comparison with Tollmien’s method will be discussed 
in the Appendix to Part ITI. 

It may be added that the adoption of the third plan leaves an error of the order of 
(aR)! in d; and @e, and an error of the order of (aR)~"* in }3 and ¢4. These errors can 
be reduced by including the higher approximations. In practice, this is hardly neces- 
sary. A detailed discussion of numerical accuracy will be found in the Appendix to 
Part III. : 

Having thus established the region of validity of the asymptotic solutions, we 
shall try to settle a few questions of considerable 
dispute, namely, (a) the “crossing substitution,” 

(b) the inner friction layers, and (c) the complex ‘ 
conjugate of the inviscid solution. ¢j>0 

a) The crossing substitution. From previous dis- y, Ya 
cussions, it is evident that if we pass from 
y>Re(yo) to y<Re(yo) along a path below the 
point yo, we are always in a region of the y-plane 
where the above asymptotic solutions hold, and no 
further investigation is necessary. In fact, ifc;>0  ¢=0 











(and is small) and Re(w¢ ) >0, the point yo is above Y, Yo Ye 
the real axis, and the asymptotic solutions are valid 

along the real axis of y. In the case of real c, the 

point Yo is on the real axis, and there is one point on 

the real axis where the asymptotic solutions fail to = 

be valid. In the case c;<0, and Re(w¢)>0, the <0 j POE y 
point yo is below the real axis, and the lines Yo ° 
arg }ao(y—yo)} = —77/6, 2/6 intersect the real Fic. 2. Diagram showing the rel- 


axis in two points y/, y/’ with* y<y/ <y/’ <y2. ative position of the real axis and the 
Thus, the asymptotic expressions (4.14) and (4.21) _ region of validity of the asymptotic 
represent one solution for yi:Sy<y/ and Solutions in each of the three cases. 

yy <ySy2, but not the same _ solution for 

y/ <y<vy/’. It is necessary to obtain a suitable “crossing substitution” in order to 
obtain the correct solutions for 7/6<arg { ao(y—yo) } <5/6 (i.e., in crossing the 
lines arg {ao(y—yo)} = —72/6, 7/6). For this purpose, we must obtain the asymp- 
totic expansion of the Hankel functions H/3[2(iaon*/?)/3], (7 =1, 2), proper to that 
region. The analytical expression for H{%} would then be quite different from that 
given in (5.2). Thus, in crossing the two points y/ and y/’ of the real axis, the asymp- 
totic solutions fail to be analytic. However, it is to be noted that the failure of 
the asymptotic solutions along the real axis does not exclude their use in the in- 
vestigation of the boundary value problems to be considered below, so long as we 
are concerned only with the eigen-value problem. It is only necessary that these 
solutions be valid in a connected region containing the end-points y; and ye. The cal- 
culation of the amplitude distribution of the disturbance (the eigen-function) in the 
neighborhood of the inner friction layers, however, is to be made with the regular 
solution, or we can calculate the eigen-function for y/ <y<y/’ by using a proper 


* The whole theory must be modified for extremely highly damped solutions for which ¥,<I, << YH» 
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“crossing substitution.” Since we are chiefly concerned with the eigen-value problem, 
we shall not go into further details. 

In order to make the situation still clearer, let us see what would happen if we try 
to obtain our solutions for y/ <y<y/’ by going along a path above the point yo. For 
simplicity, let us take the case of real c with w¢ >0, and consider the asymptotic 
expressions ¢! and ¢{ given by (5.3). We have (A, B being arbitrary constants) 


(0) —5/4 2 oe 

3; ~An exp {2 (oxen) i. (n > 0), 
(0) ,  )—5/4 3/ . —ri/4 al 

ds ~A!7| exp | 3(a0| | ) + 5ni/4}, (n < 0); 
(0) —5/4 3/2 

ds ~Bn exp }3(aon) ge (n > 0), 
(0) 1 »—6/4 3, J —5ri/4 ' 

os ~B\n| exp {#(a0| | ) + 5ri/4}, (n <0). 


These are obtained by taking a path below the point yo. If we had taken the other 
path, then arg (yn) =z for 7 <0, and we would have the functions oo” and ¢, which 
agree with ¢; and ¢, for 7>0, but are defined by 


~ (0) —5/4 | | . 3/2 32/4 


¢; ~Aln| exp {3(aoln|) e — 54i/4}, 
3/2 7xi/4 : 
— 5ri/4}, 


~ (0) 


—5/4 
$d: ~Bl\n| — exp {3(a0| 2] ) 
for »<0. Thus, if A and B are taken to be the same, we have 


~ (0) (0) ~ (0) .. (0) 

od; =— ids , 4 = — ids , for 1 < 0. 
Hence if we took $2 and $2 as the proper determinations, we would have to make 
the following “crossing substitution” corresponding to a passage from 7>0 to 7<0: 
Go id? ; ig’. If we note that 6 <«<¢™ both for w—c>0 and for w—c<0, 
we would also have the following equivalent change : 6 +6 +16 ; 62 —i¢ -G!. 
These may be compared with Eq. (16), p. 589 of Heisenberg’s paper. In making the 
comparison, we note his definition of the angle of w—c (p. 585), and the difference 
of notation in the fundamental equation of stability. 

The first study of “the crossing substitution” seems to be due to Stokes in con- 
nection with the asymptotic expansions of Bessel functions. It may therefore be 
properly designated as Stokes’ phenomenon [76]. We should also compare our re- 
sults with the work of Jeffreys [17], the W-K-B method [23] in quantum mechan- 
ics,* and the mathematical investigations of Langer [24] and others.** In those cases, 
a differential equation of the form é@’’+g¢(y)¢=0 is considered. If this equation 
is treated by the method of 4 by writing ¢ =x (n)+ex')(n) + ---,¥—yo= en, and 
q(y) =e (en) +4q0' (en)? + - , the equation for x(n) is x’’ +q0 nx =0 as com- 
pared with (4.7), xi” —inw¢ 10"? = 0, It is evident that our 7 corresponds to 77 in 
their case. Kramers has shown that the cuts in their asymptotic expansions are the 
lines arg (n) = +2/3. Thus, in our case, the cuts should be arg (n) =2/6, 52/6. This 
agrees with our previous discussions. An important difference is the following. In 


* I am indebted to Professor P. S. Epstein for calling my attention to this comparison. 
** For example, S. Goldstein, Proc. Lond. Math. Soc. (2) 28, 81-90 (1928); C. C. Hurd, Téhoku 


Math. Journ. 45, 58-68 (1939), and the papers of W. J. Trjitzinsky and others quoted there. 
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their case, the two boundary points on the real axis are separated into two regions of 
the complex plane by the cuts, so that a crossing substitution is absolutely necessary. 
In our case, the two boundary points on the real axis belong to the same region, and a 
crossing substitution is superfluous, so far as the eigen-value problem is concerned. 

b) The inner friction layers. There is also a very significant physical interpretation 
associated with the “crossing substitution” of the asymptotic solutions. The initial 
approximations ¢!” and ¢{” satisfy the inviscid equation. Hence, if c;>0, these solu- 
tions hold throughout the interval (31, ye) of the real axis, and the effect of viscosity 
is entirely negligible inside the fluid for sufficiently large Reynolds numbers. If c; $0, 
the inviscid solution can never hold all along the real axis, and hence the effect of 
viscosity inside the fluid is not negligible, however large the Reynolds number may 
be. The singularity of the asymptotic solutions means a very rapid change of velocity 
within a small distance so that the effect of viscosity is no longer negligible there. 
Physically, such a point on the real axis corresponds to a layer of fluid where the 
viscous forces play an important role. 

Referring to the foregoing discussions, we see that there are two inner friction layers 
for the damped oscillations, one for the neutral oscillations, and none for the self-excited 
oscillations. 

In the neutral case, the first term of (4.1) disappears at the critical layer w=c. 
The equation then represents a balancing of the vorticity transferred by the disturb- 
ance and that diffused away by the effect of viscosity. It is therefore understandable 
that the effect of viscosity must be predominant there. In the other two cases, w—c 
never vanishes in the fluid, there is the vorticity carried by the main flow (relative 
to an observer moving with the phase velocity c,) and there is always the change of 
vorticity due to amplification or damping. In the case of amplified oscillations, these 
two effects can be in equilibrium with the transfer of vorticity due to the disturbance, 
and the effect of viscosity is completely negligible at very large Reynolds numbers. 
In the case of damped oscillations, these effects presumably never balance each other, 
thus resulting in the formation of two critical layers, where the effect of viscosity is 
not negligible. 

c) The complex conjugate of the solution $(y). It is often argued,* that if @(y) is a 
solution of the inviscid equation with an eigen-value c, then ¢(y) is another solution 
with the eigen-value ¢, satisfying the same real boundary conditions on the real axis. 
Thus, to each damped solution, there is always a corresponding amplified solution, 
and vice versa. This argument is in direct contradiction to the foregoing discussions, 
because an amplified solution and a damped solution have entirely different charac- 
teristics with respect to inner friction layers. It appears, therefore, that ¢(y) should 
still represent a solution of the same nature as $(y). 

This is indeed the case, and can be seen more clearly from an examination of the 
complete equation (4.1). If we take the complex conjugate of that equation, and write 
y for § (which is essentially done in the usual argument), we have 


{ w(y) — c} { 6” — a6} — w"(y)b = — { d* — 2a*6” + at}. (5.5) 
a 


* Heisenberg, loc. cit. p. 596; Tollmien, loc. cit., 1935, p. 88. The failure of such an argument would 
indicate that Heisenberg’s classification of velocity profiles on p. 597 of his paper is untenable. 
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The complete stream function ’(x, y, t) satisfying this equation is 
p' (x, y, t) = B(y)e~ta (zee), 


Thus, if Im(c) <0, then Im(é)>0, and we still have a damped solution. This should 
also hold for the inviscid equation, since it is regarded as a limiting case of the viscous 
equation. From the inviscid equation itself, there is no way of telling whether 
Im(é)>0 corresponds to damping or to amplification. In fact, the asymptotic solu- 
tions of Eq. (5.5) (including the limiting inviscid solutions) hold for 


— r/6 < arg {&(y — yo)} < 74/6, w(yo) = é. (5.6) 


Thus, we have a solution ¢(y), valid in a region which is quite improper for an asymp- 
totic solution of (4.1). [Compare (5.4) and (5.6).] Hence, it is not legitimate to con- 
clude that a solution of a different nature can be obtained by taking the complex 
conjugate. The influence of these discussions upon the usual conclusions regarding 
stability in an inviscid fluid will be discussed fully in the next part of the paper. 

6. The boundary value problems. Having fully investigated the solutions of the 
equation of disturbance, we shall now turn to a 
study of the boundary-value problems which have 
been taken up briefly at the end of §3. In general, 
the boundary conditions are essentially that the 
velocities of the disturbance should vanish on the 
solid boundaries, and also at infinity if the field of 
flow extends to infinity. However, it is often con- 
venient to use equivalent boundary conditions for 
certain types of velocity distributions. (2) 

In order not to be lost in too much generality, 
we shall limit ourselves to three classes of velocity 
distributions (as specified below and shown in 
Fig. 3), and select our fundamental interval (1, ye) 
so that w’(y) 20 for v1 Sy Sye. We shall define our 
characteristic length so that ye—y.=1, and let 
oily; c, a, aR), do(y; c, a, aR), os(y; c, a, aR), 
$a(y; ¢, a, @R) be a fundamental system of solutions Fy. 3. The three types of velocity 
(4.1) arrdnged in the order discussed above. distributions. 

Case (1). Flow between solid walls in relative 
motion. In this case, the boundary conditions are given by 


wy) 


(1) 














(3) 








Pr Pp ————-——}--4— 


y, 


$(¥1) = $'(y1) = O(¥2) = $'(y2) = 0, (6.1) 


because the velocity of the disturbance should vanish on both the solid boundaries. 
The determinantal equation corresponding to these conditions is 

| dir Dor G31 Pat 
| dio boo p32 dae 
| / , , , 
| dir Gar Par ar 
} , , ’ , 
| dio hoe P32 ae 


where on, $11, etc., stand for ¢1(y1), $7 (v1), etc. In this and all later discussions, a 


| 
F(a, c, aR) = | = 0, (6.2) 
} 


-_ 
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subscript 1 or 2 attached to a function of y shall denote the value of that function at 
y=¥1 Or y= ye respectively. 

CASE (2). Symmetrical flow between solid walls at rest. In this case, it is easily seen 
from (4.1) that the disturbance can be separated into two independent parts, one 
symmetrical with respect to the line y= y2 and the other antisymmetrical. (a) If @(y) 
is a symmetrical function (antisymmetrical disturbance), then the conditions 


o(¥1) = o'(91) = O'(92) = 0" (42) = O (6.3) 
hold, and we have the determinantal equation 


$i: G21 31 ai 
$11 $21 $31 P41 





F.i(a,c,aR) =| , : ; , |=. (6.4) 
| diz Gon 2 pao | 





| is do2 32 ae 


(b) If é(y) is an odd function of y— ye (symmetrical disturbance), then the boundary 
conditions are 


o(v1) = o’(v1) = o(¥2) = 0’ (ye) = 0, (6.5) 
and we have the relation 

Gir G21 G31 a1 | 

Pi2 G22 G32 Par | 
a. ge gt ae (6.6) 
gir O21 G31 Par | 

diz 22 32 baz | 
Case (3). Flow of the boundary-layer type. In this case, the point ye is taken to 


correspond to the “edge” of the boundary layer, beyond which the velocity is sub- 
stantially constant. The boundary conditions to be satisfied at y; are the usual ones; 


F;(a, c, aR) = 


(v1) = $'(y1) = 0. (6.7) 


The boundary conditions for y becoming infinite are to be replaced as follows.* Since 
the particular integral ¢, becomes infinite as y becomes infinite, our boundary con- 
dition requires that @ is a linear combination of ¢1, ¢2, ¢3 alone. Thus, 


d = Cidi + Cope + Caos, (6.8) 


where Ci, C2, C3 are constants of integration. Also, the integral ¢; makes practically 
no contribution for y2ye so that we expect ¢$(y) to satisfy the inviscid equation for 
y2vyo. Here w’’=0, and hence two particular integrals are e+*”. The condition that 
¢—0 as y~ excludes the integral et”. Hence, ¢ must be proportional to e~*” for 
y> ye. This may be expressed as follows: 


¢ +a¢=0 for y2 yo. (6.9) 


Hence, we have the determinantal equation 


* Cf. Tietjens [72] and Tollmien [73]. 
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$11 p21 31 | 
, , | 

Fi(a, c, aR) = | bdi2 + adie $22 + Pax O = 0. (6.10) 
| ou 21 $31 | 


We note that the point ye can be replaced by any value of y>zye2. This is equivalent 
to the fact that the thickness of the boundary layer cannot be definitely defined. The 
larger this thickness is taken, the more accurate the results should be .* 
The functions F;, Fo, F;, and F; are entire functions of the parameters a, c, and R. 
Reduction of the equations for large values of aR. The equations (6.2), (6.4), and 
(6.6) can be substantially simplified for large values of aR. Referring to (4.21), we 
see that $3(y)=A(y)e~’, da(y) =B(y)e”, where A(y) and B(y) are of the order of 


unity, and FY is defined by the relation 


Y= J ViaR(w — c) dy. 


yl 


Hence, we have the following relations, giving the order of magnitude of certain 


quantities: 
$3 pon See A; tn. he ) 
ee VV iaR(w; — c) + ah one SE a eP, | 
P31 A, 031 Ay 
$32 —— 1, As | 
a = hoe. aR(we — c) ood Opin, P 
931 \ A, if | 
” f (6.11) 
32 ? Ag — 
— = ~jaR(w. — c) —+ O(WaR) pe’, 
031 t A, 
Ca | 2 
> ¥_ [iaR(we — c)}*? — + O(aR) be? 
os: is f 
a1 halaman By 42 B, 
— = ViaR(wi —c) +—) —o—e, | 
pai B, P41 B, | 
$42 a? "ae | 
— = <V/iaR(w. — c) —+— pe’, 
P41 \ oy B, 
oy : p (6.12) 
? 2 Rs a | 
== J iaR(ws —c) —+ ON aR) e”, 
oi By 
due fy, ae ae 
—— = [iaR(we — c)}*/*— + O(aR) ¢ e?, 
P41 \ By 
where 
P= f ViaR(w — c) dy. 


yl 
It then appears that the sign of the real part of P is of consequence. It can be veri- 


* In later calculation of the Blasius case, we shall take a boundary layer about 1.19 times as thick as 


that used by Tollmien. 


— 
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fied that it is always positive when c;>0. For then the path of integration can 
be taken along the real axis of y, and we have —a<arg (w—c) <0; consequently, 
—w/4<arg (P)<7/4. With reference to (4.21), (6.11) and (6.12), we see that the 
condition P = n7i (n=an integer), expresses the fact that @3i d47 =@4i 32 when terms 
of the order (aR)? are neglected. This is the corrected form of the first solution of 
Heisenberg as expressed by the condition (27) on p. 596 of his paper. Heisenberg 
also stated that such a condition can only be satisfied for damped solutions. In fact, 
from the condition just obtained for c;>0, we see that Re(P) can be negative only 
for highly damped solutions, for which the whole discussion must be modified. (Cf. 
footnote on p. 29, §5.) 

Neglecting quantities of the order e~? and (aR)— against quantities of the order 
of unity, we have the following simplifications of Eqs. (6.2), (6.4), and (6.6) for 
Cases (1) and (2). , 

CASE (1). Flow between solid walls in relative motion. We have 


Sila, ¢) es oar Gaz hola, ¢) ; 
Ss(a, c) $31 42 fala, Cc) 


CASE (2a). Antisymmetrical disturbance in a symmetrical flow between solid walls. 
We have 








(6.13) 


fola, c)/fsla, c) = 31/31. (6.14) 


CASE (2b). Symmetrical disturbance in a symmetrical flow between solid walls. We 
have 
fila, c)/fs(a, c) => ben/ bes: (6. 15) 


In these equations, the functions f;(a, c), fe(a, c), f(a, c) and f,(a@, c) are defined as 


follows: 


2 | dis | 
Pye Te tela EE Ye eed | 
| Gar ee | ber G22 | 
| ea (6. 16) 
fen OO . ide et 
| dor Poe | dor doe | 


These functions depend only on @ and c, because we may take the inviscid solutions 
for @: and ¢@2, which are accurate up to the order of (aR)—'. It may be reiterated that 
in computing 12, b22, O13, G2 in (6.16), we must take the path leading from y; to ys 
in the lower half plane. 
CasE (3). Flow in the boundary layer along a flat plate. In this case we can reduce 
(6.10) to 
fot afi ae $31 
fatafs 51 


if we also replace ¢: and ¢2 by their inviscid approximations. The equations (6.13), 
(6.14), (6.15), and (6.17) are the final equations based on which the stability in- 


(6.17) 


vestigations are to be made. 
The inviscid case. In the limit aR->«. The above equations reduce to 
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fi(a,c) = 0 for Case (1) and Case (26), (6.18) 
fo(a,c) = 0 for Case (2a), (6.19) 
fo+afi=0 for Case (3). (6.20) 


Mathematically, these are equivalent to the solution of the inviscid equation 

(w — c)(p’’.— a*d) — wo = 0 (6.21) 
subjected to one of the following three sets of boundary conditions 
$(¥1) = (v2) = 0, (91) = O'(¥2) = 0, = (91) = O'(¥2) + aG(y2) = O. (6.22) 


We have thus arrived at the conclusion that some asymptotic behavior of the stability 
conditions can be obtained by neglecting the effect of viscosity (provided proper care is 
given to the inner friction layer). This was tacitly assumed in the work of Rayleigh 
and others, while Heisenberg pointed out [loc. cit., p. 583] that a proof was neces- 
sary in accordance with some remarks of Oseen (38); he also virtually gave the proof. 

In the next part of the paper, we shall therefore consider the simpler problem of 
an inviscid fluid. After a thorough investigation of that problem, we shall investigate 
the effect of viscosity by considering Eqs. (6.13), (6.14), (6.15), and (6.17) in greater 
detail. These results may be compared with the earlier ones of Heisenberg and Toll- 
mien. 

Numerical calculations of the stability limit based upon these equations will also 
be carried out for certain important special cases. For all these purposes, the evalua- 
tion of the six functions f(a, c), fo(a, c), f3(a@, c), falar, €), 31/631, Paz/Hgg IS NECeSsary. 
We shall discuss this briefly here. 

1) Evaluation of fi(a, c), fola, c), fs(a, c), fala, c). These quantities are related to 
the inviscid solutions given by (4.14) with the path of integration subjected to the 
condition (5.4). Hence, we have 


/ / 
gu = —¢, go = 0, gil = 4, on = — —) (6.23) 


dis = (1 — c) >” a?"Ho,(c), g22 = (1 — ¢) pd a?" Kon41(C), 


n=0 

bio = (1 — oc)! DS a?"Hop_i(c) + (1 — c)— wd dre, (6.24) 

n=0 

29 = (1 — ¢)7! > a?"Kon(c) + (1 — c)~ we doe, 
n=0 
where 
Hon(c) = han(¥2), Hon1(c) = (1 — c)-*han(y2), (6.25) 
(0.29 
Konsi(¢) = Rensi(y2), Ken(c) = (1 — 6) *kans1( 92) } 


are functions of c alone. In the above evaluations, we have put w;=0, in accordance 


1945) STABILITY OF PARALLEL FLOWS 139 


with the actual conditions in all the cases considered. We have also chosen the char- 
acteristic velocity so that w.=1. Referring to (6.23), we have 


fila, c) a ser Ch22, P Sola, c) = <= Chee, 


6.26 
I3(a, c) ( ») 


, 1 ro a 
Wide2 + — diz, Sala, c) = widos + — or. 
C c 


The actual evaluation therefore depends upon the calculation of the integrals (6.25). 
2) Evaluation of b3:/63:, 42/42. These quantities are related to the highly oscil- 
lating integrals }s and ¢,4. For values of aR so large that both (a@R)"/*c, (aR)'/*(1—c) 

> 1, the approximate values of $3; /¢3, and $4./@j2 are given by using (4.21). Thus, 
$31 erties 42 gone 


(6.27) 





$5, V/aRe $1 VaR(i — c) 
The condition (a@R)'/*(1—c)>1 is generally satisfied, because c is usually small and 
aR is usually very large. This is especially.true for Eq. (6.13), for it will be seen from 
considerations of an inviscid fluid that this case is relatively stable. For the condition 
(aR)'%c> 1, the situation is different, because c is usually small. It is then more con- 
venient to approximate ¢3 by x!() given by (4.9) (or with higher approximations, 
if so desired). We then have 


31/31 = (v1 — yo)F(z), (6.28) 


where 


—2) pf i 1/2 a) 3/2 
fscod See dSH13[3(it) J 

F(z) = ————— - , 2= — ao(m — ; 6.29 
— o fre dee [3(ig)*/?] sails — 





This function has been calculated numerically by Tietjens for real values of z. We 
have made slightly more extensive and more accurate calculations. The present result 
differs slightly from that of Tietjens and is here tabulated in Table I and plotted in 
Fig. 4 together with the related function F(z) defined by 


F(z) = [1 — F(z)}-. (6.30) 
The asymptotic form of F(z) is 
F(z) ~ g73/2ertl4, (z > 1). (6.31) 


This agrees with (6.27) if a8(¥:—yo) =w¢d (v1 —yo) is replaced by ~c. 
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Fic. 4, The function #(z) shown in its real and imaginary parts (cf. Table I) 


TABLE 1.—The functions F(z) and F(z). 


| 
| 








s F, F; 7 | Fi 

1.0 0.89161 —0.35025 0.80630 —2.60557 
1.2 0.78969 —0.27310 1.77012 —2.29854 
1.4 0.71970 —0.21213 2.26836 —1.71669 
1.6 0.66931 —0.16009 2.44985 |  —1.18600 
1.8 0.63143 —0.11274 2.48104 —0.75892 
2.0 0.60144 —0.06741 2.43927 —0.41253 
23 0.57599 —0.02226 2.35196 | 0.12348 
2.4 0.55230 —0.02395 2.22724 0.11916 
2.6 0.52773 —0.07203 2.06929 0.31558 
2.8 0.49952 0.12220 1.88566 0.46043 
3.0 0.46456 0.17391 1.68938 | 0.54872 
3.2 0.41947 0.22520 1.49726 | 0.58082 
3.4 0.36110 0.27193 1.32516 0.56401 
3.6 0.28802 0.30705 1.18429 0.51074 
3.8 0.20352 0.32130 1.07982 0.43560 
4.0 0.11800 0.30721 1.01118 0.35220 
4.2 0.04698 0.26559 0.97361 | 0.27133 
4.4 0.00240 0.20811 0.96056 0.20038 
4.6 0.02160 0.14475 0.95989 0.13601 
4.8 0.01477 0 0.97659 0.09503 











(To be continued) 
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METHODS OF REPRESENTING THE PROPERTIES 
OF VISCOELASTIC MATERIALS* 


BY 
T. ALFREY 


Monsanto Chemical Company, Research Department, Springfield, Mass. 


Introduction. In a recent paper! it has been shown that the solution of the first 
and second boundary value problem for linear viscoelastic media can be obtained in 
two steps requiring (a) the solution of an equivalent problem for a perfectly elastic 
medium, and (b) the determination of the response of the viscoelastic material to an 
applied shearing stress (or shearing strain) which is a given function of time. The 
study of the behaviour of viscoelastic materials in pure shear is accordingly seen to 
be of particular importance. To coordinate various manners of describing this be- 


haviour is the purpose of the present paper. 


From the mathematical point of view the behaviour of a viscoelastic material in 
pure shear is represented by a differential relation between the shear stress s and the 


shearing strain e. We may write this relation in the form 





Ps = 20«, 
where the differential operators P and Q are defined by 
fa) m om} 
P = Papua + Pm-1 + Pree + Po, 
at ot™ 1 
fA) n o"- 1 
O= ee are + i... +-+++@. 
q oi” q oi"! 


The m+n+1 coefficients pPmt,-+:, Poy Qn***» Qo are 
constants characterizing the mechanical properties of the 
material. Equation (1) can also-be considered as the gen- 
eral stress strain relation of an incompressible viscoelastic 
medium. In this case, € may be taken as denoting any 
component of the strain tensor and s as denoting the cor- 
responding component of the deviatoric part of the stress 
tensor. 

While Eq. (1) gives a complete mathematical descrip- 
tion of the mechanical behaviour of a viscoelastic material 
in pure shear, it is often found useful to express this be- 
haviour in terms of a mechanical analogue, or model, con- 
sisting of springs and dashpots. Figures 1 and 2 show 
typical models of this kind. 

Models of the first type, shown in Fig. 1, consist of re- 
tarded elements (Voigt elements) coupled in series. Each 


* Received Nov. 15, 1944. 
1 T, Alfrey, Quarterly of Appl. Math. 2, 113-119 (1944). 
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Fig. 1. Mechanical model: 
3 Voigt elements in series. 
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element is made up of a spring coupled in parallel with a dashpot. In such a model 
the total extension (corresponding to the strain €) consists of m contributions, one 
from each of the » Voigt elements. The extension €; contributed by the ith element 
is connected with the load s by means of the relation 


s = 26:6; + 2n:é;, (2) 


where G; is the spring constant and 7; the dashpot constant of the ith element, and 
the dot indicates differentiation with respect to time. The load s is the same for all 
elements coupled in series, and corresponds to the stress in the viscoelastic body. The 
mechanical behaviour of the model is defined by m equations of the form (2) in con- 
junction with the relation e=).€; which defines the resulting extension e. 
Models of the second type, shown in Fig. 2, consist of another kind of composite 
elements (Maxwell elements) coupled in parallel. Each element is made up of a spring 
coupled in series with a dashpot. In such a model the 
total load (corresponding to the stress) is divided 





===. among the m elements. The load s; carried by the ith 
element is connected with the extension ¢ by means 
of the relation 
S 

G, G, x 1 1 
i= —i + — Si, (3) 

2G; 2n; 
Te where G; and 7; have the same meaning as above. 


7 VA The extension ¢ is the same for all elements coupled 
in parallel and corresponds to the strain of the visco- 
elastic material. The mechanical behaviour of the 








? model is defined by equations of the form (3) to- 
v gether with the relation s =)s; which defines the 

Fic. 2. Mechanical model: resulting load s. 
3 Maxwell elements in parallel. In a study of molecular mechanisms of visco- 


elastic deformation, a model of the type shown in 
Fig. 1 may be preferable to the general stress-strain relation (1). In such a study, 
each contribution to the strain may often be identified with some specific molecular 
process, and hence the strain contributions €, €, €3,° °°, €,, aS well as the total 
strain €, can be said to possess a physical significance. Likewise some authors have 
attempted to identify the various stress contributions of a model of the type shown 
in Fig. 2 with individual “molecular mechanisms of supporting stress.” From the 
point of view of mechanics of continua, on the other hand, the formulation (1) is 
preferable to any mechanical model, since in any macroscopic study only the total 
stress and total strain are observable quantities. 

If the molecular and the macroscopic methods of approach to viscoelastic be- 
haviour are not to become isolated from one another, it must be possible to change 
readily from one method of description to the other. It is the purpose of this paper 
to provide simple techniques for these conversions. The paper is divided into four 
parts corresponding to the following problems: 

1. Given the constants occurring in the stress-strain relation (1), to compute the 


constants of the equivalent Voigt model. 
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2. Given the constants occurring in the stress-strain relation (1), to compute the 
constants of the equivalent Maxwell model. 

3. Given the constants of a Voigt model, to compute the constants of the equiva- 
lent stress-strain relation. 

4. Given the constants of a Maxwell model, to compute the constants of the 
equivalent stress-strain relation. 

1. Determination of the constants of the Voigt model. 

A. Nondegenerate case. In the standard or nondegenerate form of the stress-strain 
relation (1), the operator P is of an order one less than that of Q. The relation (1) thus 
has the form 

a”—'s os 0”e . 
par + pn—2 aot + + pos = 2qn - + + 2qoe. (4) 





If both the coefficients go and q, do not vanish, the corresponding mechanical model 
will consist of m Voigt elements, all nondegenerate. If g,=0, one element of the model 
consists of a spring only, and if gy=0 one element consists of a dashpot only. These 
degenerate cases will be considered in the following sections. Cases are also possible 
where some other coefficient vanishes. This does not affect the form of the resulting 
model or the nature of the mathematical treatment. 

A given Voigt element is defined by its constants G and n. The compliance J is 
defined as the reciprocal of G; J=1/G. The retardation time 7 of the element is de- 
fined as r=»/G=Jn. Our problem is to compute, from the 2” coefficients of the non- 
degenerate stress-strain relation the 2m parameters of the mechanical model. The 
method given below depends upon the fact that both the model and the stress-strain 
relation must give the same prediction as to how the total strain will change with 
time when a given stress s(t) is applied. It is sufficient to equate the responses to the 
particular stress s(t)=t"-'. The general solution of the equation P(t"—') =2Qe is 
the sum of the general solution of the associated homogeneous equation Qe = 0 and the 
particular polynomial solution of the complete equation. In the same way, the re- 
sponse of the model to the stress ¢"~' is the sum of the general response to a zero stress 
and a particular polynomial response to the stress t"—'. If the response of the model 
is to be identical with that predicted by the stress-strain relation, the constants 
of the model must satisfy certain conditions. First, the retardation times 17; 


(i=1, 2,---,m) of the Voigt elements are the negative reciprocals of the roots x; 
of the characteristic equation g,x"+qnix""'+ - - - +Q0x+qo=0; 
1 all 
Ce -— 6 (5) 
Xs 


Thus, the » retardation times of the model are determined by the general solution 
of the homogeneous differential equation. 

In order to complete the specification of the model the particular polynomial solu- 
tion must now be used. The particular polynomial solution of the equation P(¢"~') 


= 20¢ will be of the form 
2e(t) = ao + ayt + aol? +--+ + ay". (6) 


The coefficients ao, a1, - - - , @,-1 are determined in the usual manner. 
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In order to obtain the model strain €¢ corresponding to the stress s =¢"~', we con- 
sider first the behaviour of a single element under this stress. Equation (2) can be 
written in the form 

J is = Qe; + 27;6;. 
Setting s="~—' and determining the polynomial solution of this differential equation 


for €;, we find 


ee eee ee oe per. (7) 


Since the total strain is =) €;, comparison of (6) and (7) shows that the compliances 
J; must satisfy the linear equations 


n 
oy = be J is 


t=1 


ll 


An—2 (n — 1)>0 J; Xi, 


> 2 (8) 
Qn—-3 = (n — 1)(n — 2) Ii Xiy 
n a—1 
ao = (n — 1)! >> J ;/x; ’ 


| 
| 
| 
| 
| 
| 
t=] | 
| 
| 
} 


t=1 


where the retardation times 7; have been expressed in terms of the roots x; of the char- 
acteristic equation in accordance with (5). 

The Voigt model is completely specified when the compliance and retardation 
time of each Voigt element are determined. 

B. Degenerate case; qo=0. If the coefficient go of the differential operator Q is zero, 
one of the roots, x; say, of the characteristic equation vanishes. This indicates that the 
spring constant G of the first element is zero, the element consisting of a dashpot only. 
In this case the compliances J; can no longer be found from the linear equations (8) 
because of the infinite terms 1/x,, 1/x7, - - -. A parameter may be substituted for 
the zero coefficient go, the 2” parameters 7, Gi, 72, G2, - - - , Mn, Gn, may be determined 
in terms of this parameter, and finally the parameter may be allowed to approach 
zero and the limiting values of these 2” parameters obtained. However, this involves 
a complicated procedure even in simple cases. The following alternative treatment 
of this degenerated case seems preferable. 

The total extension e of the model consists of the extension €, of the degenerate 
first element and the extension e’ of the chain of »—1 nondegenerate elements; 
€=€,+e’. The mechanical behaviour of the degenerate element is given by 


— jj = 2€1, (9) 
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and that of the chain of m—1 nondegenerate elements, by a relation of the form 
P’s = 20’, (10) 
where the differential operators P’ and Q’ are of the orders n—2 and n—1, respec- 


tively. Applying the operator Q’ to both sides of (9), differentiating (10) with respect 
to time and adding, we obtain 








1 
P's + —Q's = 20'( + &) = 20%. (11) 
1 
With 
pe o”™ 2 4 , on-3 “ ‘“ , ) 
~ op Poe - | 
r (12) 
; : o” 1 P n—2 , 
QO = qr-1 oe + Qn-2 or >>: $i. 


qn. , an 1 : Qn—2 an 2 ; qi F;) qo 
1+ “~~ )— + (pt. + “= )—— + --- + (6 +-)—+=]s 
m1 oi" m 7 ot"~* m/Obo om 
o” o=! 0 
== 2 E J + Qn-2 a + eee + qo < Ie (13) 


ot” an" Ot 


When both sides of (13) are divided by the coefficient of the highest order term on the 
left-hand side, this equation must be identical with the stress-strain relation (1) in 
which go=0. Comparison of the lowest order terms leads to the relation 


m = 9:/ po; 
comparison of the highest order terms leads to 
’ me | 
Qn-1 po 
i+—= E —-—@Qnl . 
m1 q1 


Abbreviating this expression by r, we find by further comparison of coefficients in (13) 
and (1) that 


, 


, = 
Qn—1 = 74n, In—2 = 1On—-1,°°*° 1 Go = 701 


po 
Pn-3 = r( pes ee Qn-1 ’ 
q1 


pels r(p, aes... tn-2), (14) 
q1 


Po 
r(p, a. an). 
q1 


The differential operators (12) are thus determined, and the procedure outlined under 
1A permits the determination of the constants of the »—1 nondegenerate elements. 





po 








148 T. ALFREY [Vol. III, No. 2 


C. Degenerate case; g,=0. If the coefficient g, of the operator Q is zero, one Voigt 
element of the corresponding model consists of a spring only. The compliance of this 
isolated spring can be shown to equal 1/g,-1. A procedure similar to the one developed 
above will permit to determine the constants of the nondegenerate elements. 

2. Determination of the constants of the Maxwell model. 

A. Nondegenerated case; A nondegenerate model of the Maxwell type corresponds 
to stress-strain relation (1) in which g,=0 and go=0 (i.e., to a doubly degenerate 
Voigt model). When the operators are of this standard form, the model will consist 
of m Maxwell elements in parallel. The 2m constants of this model can be computed 
from the 2m coefficients of the stress-strain relation by a method which, except for 
the interchange of stress and strain, is almost identical with that of Section 1. 

For any given imposed strain sequence e(t), the stress must vary in a definite 
fashion s(t). The predictions of Eq. (1) and the set of differential equations (3) must 
be identical for every case—in particular, for the strain sequence e(t) =¢". The results 
are as follows: 

The m relaxation times of the m Maxwell elements are the negative reciprecals of 


the m roots of the characteristic equation x"+pnix" '+ --- +pix+po=0; 
1 — 
i= - (15) 
Xi 


The specification of the model is completed by determination of the m dashpot con- 


stants 71: - - Mm. These are obtained by solving the following set of m linear equations. 


” 

a@n-1 = m >. Nis | 
ae | 

. | 
m | 

| 


In» = m(m — 1)>> ni/ Xis 


i=1 


(16) 


m 


, m—l 
ao = m! z. ie ae 
i=l 


where the a; are the coefficients of the particular polynomial solution 
pi 


s(t) = do + ayt + aol? + +++ + Anil” 


of the differential equation 
Ps = 20%". 

B. Degenerate case; 9,0, qo=0. If the order of the operator Q is one greater than 
that of P (i.e., if ¢,0), then one Maxwell element of the model consists of a dashpot 
only. The constant of this isolated dashpot is found to equal g,. A procedure patterned 
on that of Section 1B will permit the determination of the constants of the remaining 
nondegenerate elements. 

C. Degenerate case; qo#0, gn=0. If the coefficient go does not vanish, then the 
model contains one element which consists of a spring only. The constant G, of this 
spring, is found to equal go/po. The constants of the remaining nondegenerate ele- 
ments can again be determined by a procedure similar to that of Section 1B. 
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3. Determination of the operators P and Q from the constants of a Voigt model. 
Consider a material whose behavior in shear is reproduced by a model consisting 
of m Voigt elements in series. The 2” parameters of this model are known. The equiva- 
lent relationship between stress and strain can be determined by either of two 


straightforward methods. 
1. The method of part 1A can be used in reverse. This immediately gives an opera- 


tor which is directly proportional to Q. 


"(a 1 
dQ = II(=+ ~), (17) 


i=] Ti 
where A is an undetermined multiplier. 
The operator P can subsequently be determined by equating the particular poly- 


nomial responses to a stress s=?t"—'. 
2. The mechanical behaviour of the Voigt model is expressed by the following set 


of equations: 
s = 2Gye, + 2m16€1, 


s = 2Gre2 + 2n2€2, 


; ; > (18 
s5= 2Gn€n + 2nn€ns 


n 
e=) €;. | 


i=1 





The uth equation can be rewritten, as 


n—1 n—1 
s= 2G, (. -> «) + 2nn (« -> a), (19) 
i=] i=] 

If each of these equation is differentiated (n—1) times, a total of m? equations will 
result. These equations will contain (m*—1) derivatives of the form 0’e;/dt’. All of 
these derivatives can be eliminated, leaving a differential relation between s and e, 
by multiplying each of the n? equations by an appropriate factor and adding. The 
determination of the factors may, of course, be rather cumbersome. 

3. The problem can, however, be simplified by a judicious combination of meth- 
ods (1) and (2). We determine first the operator AQ in accordance with (17). We then 
formulate the set of m? equations considered above. m of the necessary n* factors can 
immediately be written down. They are obtained from the coefficients of the opera- 
tor (17). The form of the n? equations is such that the remaining factors can be evalu- 
ated one at a time if the above set of m factors is known. The result of this procedure 
is the desired operator equation. 

4. Determination of the operators P and Q from the constants of a Maxwell 
Model. Consider a material whose behaviour in shear can be reproduced by a model 
consisting of » Maxwell elements in parallel. The 2” parameters of this model are 
known. The equivalent relationship Ps =2Qe can be determined by methods almost 
identical with those of Section 3. Only the simplified third method will be repeated 
here. 

The operator P is given by the equation 
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. 0 1 
AP = |] (; + -), (20) 
imi \OL Ti 
where X is again an undetermined multiplier. The mechanical behaviour of the model 
is expressed by the equations 
1 1 ) 


é=—h+—s: 
2G, 2m) 


aie Rha ; - 
°" 2G, (: ~ 2 s:) . 2nn (: 7 X s) 


If each of these equations is differentiated (n—1) times, a total of m? equations are 
obtained, involving n?—1 derivatives of the form 0’s/d?’. All of these derivatives can 
be eliminated, leaving the desired stress-strain relation, by multiplying each equation 


(21) 


by an appropriate factor and adding. The » coefficients of the operator (20) provide 
of these factors. The remaining (n?—™m) factors can then be obtained one at a time. 


151 


RAYLEIGH WAVES AND FREE SURFACE REFLECTIONS* 


BY 


C. H. DIX, C. Y. FU, ann ETHEL W. McLEMORE 
United Geophysical Company, Pasadena, California 


1. Introduction. The theory of waves associated with the plane boundary of a 
semi-infinite, isotropic, homogeneous, perfectly elastic medium was first given by 
Lord Rayleigh,’ who discussed the problem for plane waves of fixed frequency. Many 
papers have been written giving treatments of variations of the problem studied by 
Rayleigh but the treatment in Rayleigh’s original paper contained most of the results 
of interest for plane waves. 

Had Lord Rayleigh realized the great practical importance of his surface waves, 
he would doubtless have included more numerical results in his original paper, and 
the material of the present paper would have been more or less completely included 
therein. Rayleigh waves are important in the seismic method of oil exploration 
since they generally occur as a troublesome noise on reflection seismograms. 

In the present paper, we are interested in the theory of the reflection of a plane 
compressional incident wave at the free surface. In this theory a cubic expression 
occurs which also occurs in the theory of the Rayleigh waves. Furthermore our inter- 
est is primarily in obtaining numerical results, so that examples may be readily pic- 
tured. The results are primarily of theoretical interest since our waves are never plane 
and the medium is only rarely approximately homogeneous. 

2. Reflection at the free surface. This problem originally treated by Knott? and 
Zoeppritz leads to the relation 

R v* sin 27; sin 2i — V? cos* 2r; 


, 
I v? sin 2r,; sin 27 + V? cos? 2r; 





where R, J, m, i, v, and V are respectively the amplitude of the reflected compres- 
sional wave, the amplitude of the incident compressional wave, the reflection angle 
of the shear wave, the angle of incidence, the velocity of the shear wave, and the ve- 
locity of the compressional wave. 

We may make the substitutions w=sin? r;= p*v? and s=A/u=20/(1—2¢) and ob- 
tain . 


R_— 4w(1 — w)'/2[1 — (s + 2)w]"/* — (s + 2)%(1 — 2w)? 





I 4w(1 — w)"2[1 — (s + 2)w]*/? + (s + 2)"/2(1 — 2w)? 
N(s,w) N? ND 
D2 








D(s, w) ND 


(1) 


Now real values of both N(s, w) and D(s, w) are graphed on Fig. 1, for various values 


* Received Dec. 2, 1944. 

1 Lord Rayleigh, On waves propagated along the plane surface of an elastic solid, Proc. Lond. Math. 
Soc. 17, 4-11 (1887). 

2C. G. Knott, Reflection and refraction of elastic waves, Phil. Mag. (5) 48, 64-97 (1899). 

3K. Zoeppritz, Uber Erdbebenwellen VIIB, Géttinger Nachrichten 1919, 66-84. 
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of s or a. It is seen that the graphs break up into two sets, one for w31/(s+2) and 
the other for w21. When w=1/(s+2) or w=1, the graphs of N and D both have 
vertical tangents. It will be observed that N has two or no real zeros between 0 
and 1/(s+2) while D has one real zero 21.09574. The double zero corresponds to 
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™ Fic. 1. 


o =0.26308207, w =0.27969015. For large values of wwe have D~ —2(s+1)(s+2)7'/*w 
and N~ —8(s+2)'/2w?. Hence the real zeros for w>0 are fully accounted for. 

There is some interest in the plot of R/J against i for various a’s. This is shown in 
Fig. 2. Observe that the tangent to the curves at 1=0° is horizontal and is vertical 
at i=90°. Observe also that a discontinuity exists for o=0 at i=90° since 
lim;.90°(R/I) emo = +1 and (R/1)i<90°,0<e<1/2= —1. One may note that the zeros of 
(R/I) correspond to i’s, for the given o, for which there is no reflection of energy in 
the wave of compressional type. The change of sign of (R/J) corresponds to a phase 


reversal. 
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© = 26306207 














Fic. 2. 


Since the zeros of N(s, w) are of considerable interest we have plotted io against o 
(where ip) is incidence angle for which 
N=R=0) in Fig. 3. This graph brings out 





a point which is indeed curious, namely if 
o =0.15 then for 7p = 87.76° all the reflected 
energy appears in the shear wave, whereas, 
if we add only 2.24° to z all the reflected 
energy appears in the compressional wave. 
Birch* records Poisson ratios for granite 
blocks of 0.093, 0.096, 0.116, 0.086 and 
0.109. These are selected low values. It 
is established in a later paragraph that 
(dio/do),.0.= — © for the upper branch 
and (dio/do),-o=+0.178 for the lower 
branch and (dio/da) 0.23 = © for the upper 
and lower branches. 


4 Francis Birch, Handbook of physical constants, 
Special Paper no. 36, Geol, Soc. of Amer., pp. 73-74 
(1942), Fic. 3. 
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If we rationalize either the denominator or numerator of (1) as indicated we obtain 
y = ND = 16(s + 1) w* — 8(3s + 4) w? + 8(s + 2)w — (s + 2). (2) 


The curves for y=y(s, w) are plotted in Fig. 4. The fixed point at w=1.0957444, 





ND? Y 2; 6:4 


» od 








W = 1.0957444 
Y = -1.83927 














Fic, 4. 
y = — 1.83927 will be noted. It corresponds to the w for which dy/ds =16w* —24w? 
+8w—1=0. Note also that 
y(s1, w) — y(Se, w) = ($1 — S2)(Oy/9s). (3) 


If either N or D is zero then y is zero. Assume y=0 while s and w vary. Then after 
a few reductions we obtain 


dig 9 (s + 1)?(s + 2) (dy/ds) 





(4) 


a sin 2ée (8y/0w) wav, 
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This relationship should be kept in mind in connection with the plot of Fig. 2. Note 
that the double zero of the cubic (2) is the common zero of y=0 and dy/dw=0 (yield- 
ing after elimination of w the relation 33s*+12s?—27s —30 =0, so that sp = 1.11043541 
and ¢ =0.26308207 while ip = 68°51’44’’;5 hence the 0y/dw in (4) is zero, leading to 
the vertical tangent. Note that the point ¢=0, ip = 90° is not attained on Fig. 3 but 
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that as we approach this point from ¢>0, dio/da—>— © because sin 2i9—0 where- 
as the other factors are bounded away from zero. 

3. Rayleigh waves. In the usual theory of Rayleigh waves® a cubic equation oc- 
curs which, using our s-notation, can be written in the form, 


16(s + 1) — 8(3s + 4) + 8(s + 2)w? — (s + 2)0* = 0, (5) 


where #= V32/v’, v and Ve being respectively the velocities of the shear wave and 
Rayleigh wave. Thus when y=0, w=a%~-'. But w=sin? 7,=p*v? where p=sin i/V 
=sin 7,/v. Hence p*v? =v*/VRz so p=1/Ve. Thus p, which for real values of i or n1, 
can be interpreted in terms of the reciprocal of the velocity with which the wave 


5 B. Gutenberg, Energy ratio of reflected and refracted seismic waves, Bull. Seis. Soc. Amer. 34, 85-102 


(1944). 
* J. B. Macelwane, Theoretical seismology, Part I, Wiley & Sons, New York, 1936, p. 114, Eq. (5.33). 
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sweeps along the surface, 0¢/dx, has the same interpretation in the case of the Ray- 
leigh wave. 

The Rayleigh wave case corresponds to the zero of D(s, w) which in turn yields 
a poristic problem when J=0.’ 

The Rayleigh wave velocities can be readily computed by solving the cubic 
equation or by an inspection of Fig. 4. However, Figs. 5 and 6 show V, v, and 
Vr for the respective cases where 4 =constant = 2.96 X10" dynes/cm.?, and where 
k = (u/3)(3s+2) =constant = 4.9410" dynes/cm.?=modulus of compression for 
various values of o. 

7T. Sakai, On the propagation of tremors over plane surface, Geophysical Magazine, Tokyo, 8, 1-71 
(1934). 
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ON THE NON-LINEAR VIBRATION PROBLEM 
OF THE ELASTIC STRING* 


BY 


G. F. CARRIER 
Harvard University 


1. Introduction. It is well known that the classical linearized analysis of the vi- 
brating string can lead to results which are reasonably accurate only when the mini- 
mum (rest position) tension and the displacements 
are of such magnitude that the relative change in | 
tension during the motion is small. The following Ty] \ YO 
analysis of the free vibrations of the string with fixed 
ends leads to a solution of the problem which ade- Vix) 8 U(x+4x) 
quately describes those motions for which the changes 
in tension are not small. The perturbation method is 
adopted, using as a parameter a quantity which is T(x+ax) —~ 
essentially the amplitude of the motion. The periodic 
motions arising from initial sinusoidal deformations Fic. 1. Displaced element of string. 
are closely approximated in closed form. The method 
is applied to motions not restricted to a single plane and finally the exact solution for 
the transmission of a localized deformation is indicated. 

2. The equations of motion. The equations of dynamic equilibrium of an element 
of the string, deformed into a plane curve as shown in Fig. 1, are 
0*v 


é [7 cos@] = pA (1) 
’ —~is 0080) = pA =——* 
Ox ‘ Or? 


Ax 








O7u 





2 


re) 
— [T sin @] = pA 
Ox 


where p denotes the mass per unit volume, A the cross-sectional area of the string in 
the rest position, and @=arc tan [u’/(1+v’)], the primes indicating differentiation 
with respect to x. The condition of fixed ends implies that, 


l 
f v'dx = 0 for allt. (2) 
0 


The stress-strain relation of the string is assumed in the form, 

T — T) = EA{[(1 + 0)? + (w’)?}!2 — 1}, (3) 
where 7» is the tension in the rest position and E is a constant characteristic of the 
string material. The following dimensionless quantities are introduced to simplify the 
algebraic work 
To T — T> : Wx T (=) 

= ’ a=) =) 9 2 oj ——— ‘ 
EA To l l \pA 


After differentiating Eqs. (1) with respect to x, setting 


* Received Jan. 3, 1945. 
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u’ u’ 
sin 0 = — —=~y, 
[(w)?+(+v)?] itaxr ° 
cos 8 = (1 — y?)}, 
and eliminating v’ between Eqs. (1) and (2), we obtain 
g2 92 . 
a [a + De] = = [A + ae], (4a) 
2 n? 
0? , : 0? . 
ry [(1 + 7)(1 — ¢*)¥?] = On [(1 + a?r)(1 — g*)*/?], (4b) 
f (1 + a?r)(1 — g*)2%dt = x. (4c) 
0 


These equations rigorously define the motion of the string which is acted on by on 
external forces. 
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Fic. 2. Comparison of periods obtained by - 
linear and non-linear theories. Fic. 3. Period v.s. initial tension.* 
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Motion defined by Eqs. (15). a0. 


3. The perturbation procedure. It is convenient to choose, as the perturbation 
parameter of the problem, a number e€ which is essentially the amplitude of the mo- 
tion.** The two functions ¢g and 7 are therefore expanded in powers of this parameter 


as follows: 
g=aleyiteystey +--- ], tr=eretemt:::. (5) 


It is easily seen that a reversal of the sign of e should merely reverse the sign of ¢. 
Hence the omission of the even powers of ¢ is justified. In a similar manner the fun- 
tions 71, T3, -- - can be seen to vanish. That ro vanishes is seen by inspection of 
Eq. (4c). The expressions for g and 7 are now substituted into Eqs. (4), the coeffi- 


* In Fig 3. the ordinates should be labeled P/P*. 
** Equation (15) indicates more precisely the meaning of «. 
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cients of each power of € are equated to zero, and the following system of equations 
is obtained: 





2 
Lo(¢1) _ 0, (6a), Li(r2) caidas a*Lo (=), (6b) 
24 
, a 1 
Lo(¢s) _ Ly(r2¢1), (6c), Ly(r4) _ a*Lo (ose: a 8 ). (6d) 


2 


3 
. g $1? g 
Lo(¢s) a Li(tayr + T2P3), (6e), Li(r6) = — aL (ov + —t a? : : + a=), (6f) 





. 
, 


where 
0? 0? 0? 0? 
L=—-—> 1=a°-———- -——») 
0g On? On? og 
and 
: 2 
gy 
f (". - *) dt = 0, (7a) 
0 2 
2 4 
9 a $1 
f (". — gigs + ar ) ae = @, (7b) 
0 


Since each of the operators in the foregoing equations is linear, it is now a simple 
matter to evaluate successively the g; and the 7;. For the moment, we confine our 
attention to the motion defined by choosing as a solution to Eq. (6a) the function, 


¢1 = cos £ Cos 7. (8a) 
Note that for g,=cos né cos nn the same solution will exist when / is replaced by 1/n 
in the definitions of £ and 7. Solving successively Eqs. (6), starting with the foregoing 
definition of g:, and using Eqs. (7) to determine the arbitrary terms appearing in the 


T;, we obtain 





1 a? 
t?~= 7 cos? » + 7 cos 2&, (8b) 
‘as bavi ‘| - lth dn EO EG dos. FO cos | 
32 128 128 
a?(9 — a?) 
—— -— cos 3£ cos n, (8c) 
3 — 2a? — at 1 — 9a? 3a? 1 
™ = | -— = n sin 2n + —— (cos 2n + cos 4n) — cs. * - oo" r| 
a4(21 — a? 3a? 13 — a? 
—_—_—_—__—_—_—_ (08 2é +e —— oer cos 4é cos 2n, (8d) 


“512 2048 4—@ 
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71+a*+:---) 5 ; 9 ; in 

¢s = cos £| — —————___——— 7? cos n — —— n sin n — —— nsin 3n + 27" cos Sn | 
2048 512 4096 

Be aie (8e) 


The arbitrary solutions of Eq. (6a) which may be added to each of the ¢g; as they are 
evaluated have been chosen in such a manner that lim ae’g; exists when a tends to 
zero and ae=constant=a.* This limiting process, of course, defines the motion 
wherein the initial tension TJ) is zero and the amplitude a is non-vanishing. An in- 
vestigation of this problem will simplify the question of the convergence of the func- 
tions ¢ and 7 as defined by Eqs. (5) and (8). When a tends to zero as specified above, 
the symbols 7 and » become meaningless. Hence, we replace them by 


o¢=—=$a’r, and 7 = as. 


The limiting process then yields the following expressions for the g; and the g; 


aey; = acosé, 


iy~s\ 9 9 
ae’y3; = a*| — : cos § +- cos — ~ —— cos 3é |, 
¥. 128 


, Tass Ife 45 /s\? 
ae’g, = a° ( cos & + ( ) cos & — cos 3& + f(€) |, (9) 
4!\ 2 ae 128\ 2 . 


1 
€°09 =- a°*, 
4 | 
1 cos 2 37 
eto, = a*| —— 5? + —— 4 — |, | 
16 32 256 
Pes 13s? 3s? 
e6o, = a® — + — cos 2¢ 
128 512 128 + (10) 
1 ¥1 6 : 
+ er + so9i¢3 + a2 — |a*e® + g(E), 
56 7 . 
e8g3 = a®| —— + ; 
1280 


These solutions may also be obtained, of course, by assuming @ equal to zero at 





* Such complementary solutions are usually chosen to be consistent with a given set of initial condi- 
tions. However, it is convenient here to choose them so that the solution does not become meaningless 
when a0, ae=a. Equations (15) indicate that this choice leads to a solution corresponding to a nearly 
sinusoidal initial deformation. 
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the outset, expanding ¢ and ¢ in powers of a parameter a, and proceeding in the fore- 
going manner. 

Note that the leading terms of the g; define an absolutely converging series for 
all a. Note also that the remaining terms of each g; are dominated by this leading 
term. In fact, for sufficiently large s, the sum of the remaining terms in each g; is 
as small as we please compared to this leading term. Although this dominance has 
not been shown to occur uniformly, it is to be expected that the series defined by 
Eqs. (9) and (10) will converge over some range of a. The requirement, “sufficiently 
large s” introduces no difficulty since the initial value of s may be chosen arbitrarily 
large. 

The functions g and o are now most conveniently written in the forms 


g(, s, a) = afi(as, £) + a*fs(as,&) +---, 
a(t, s, a) = a°go(as, £) + atgs(as, =) +---, 
where the terms of the series defining the f; and the g; are easily chosen from Eggs. (9) 


and (10). f: and ge are composed of the previously mentioned leading terms, and it is 
easily established that they converge to the values 


(< 1 ) 1 (< 1 ) (12) 
f,=cn{—, —=)cosé, go = —cn?{—, —=}, 
" 2 V2 . i 2 V2 


where cn denotes the elliptic cosine. Energy considerations may be used to show that 
the remaining f; and g; are bounded, and it is to be expected that the motion is closely 
described by ¢=af; and o=a*g. when a is sufficiently small. For most materials, a 
value of a? greatly in excess of 10~ will lead to plastic deformations; hence, the motion 
of such strings is well defined. 

The motions arising when 7» is arbitrary, as defined by Eqs. (5) and (8), can also 


be written in the form, 
aeF )(€, Up €) + a®'eb F3(é, n; €) + et + P(é, n, a, €), 
&(G(E, n, €) + a%’Ga(é,n,e +-°> | + QE, , a, €), 


(11) 


II 


13) 


7 


where P and Q are those parts of g and rt which vanish when a tends to zero and 


ae=a. For this case, 


i og: 1 Se 
F, = cn (i+—~%, k) cos g, G. = — cn (y/1+5 #), 
1=¢ (y ri U] ( 4 4 ’ (14) 


where k =¢€[2(4+e?) | 2. It is evident, in view of the foregoing results, that 


lim F ;(é, n, €) _ SiG, as) 


and it is to be concluded that since the series defining the F; converge as € tends to 
infinity, they will also converge for the smaller values of «. Both a and ae must be 
small because of elastic considerations, which indicates that P and Q will also exist. 
We conclude therefore that the motion of the string, whose “amplitude” ae is of the 
order of magnitude required by elastic considerations, is adequately defined by the 
leading terms of Eq. (13). That is, in the first approximation, 
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aeen| 4/1+~<, & | cos 
é? é 
raion | 4/14 <n, |. 
4 4 


Figs. 2 and 3 compare the results of this analysis with those of the linear theory. 

4. The motion following an arbitrary initial deformation. The motions derived in 
the preceding section are obviously those corresponding to initial sinusoidal deforma- 
tions. If the perturbation procedure is again carried out, and if for g; the function 
¢1=)_;b; cos j€ cos jn is selected, a solution will be obtained, the leading terms of 
which contain no powers of @ greater than unity. The solution so obtained will corre- 
spond to an initial deformation, ¢;(£, 0) =>_;b; cos j#£. This predominating part of the 
solution may, however, be obtained by a simpler, less rigorous, procedure which 
nevertheless leads to identical results. We merely expand (1—g?)!/? in the conven- 
tional power series and omit in Eqs. (4), ¢"*? as compared to ¢", and a? as compared 
to 1. We thus obtain as replacement for Eqs. (4) 


6 
Il 


(15) 











0? dy 
ag? pet tel a | 
Or 
38 = 0, hence r= 7(n), | (16) 
r ¢? 
27 — — jdt = 0. 
i} (« , + ’ 
Finally the first of these becomes 
an? 7 dy 00 
1+——f 0*(E, |e =. 17 
| Te Pp ¢7(&, n)dé ae an? (17) 


The solution corresponding to the initial conditions specified at the beginning of this 
section is found by considering that solution of the form g=a)_; b; cos jéy;(n), where 


y;(0) =1 for each 7. 
Upon substitution of this function, Eq. (17) yields the following set of ordinary 


differential equations 
2 1 2.2 
Vn’ + n va[1 + 4 , A = 0. (18) 
3 
These may be written in the conventional operational form 
2 2 n* 9 
(D +n hn = — ride bi, (19) 
j 

and standard integration procedure leads immediately to the integral equation 


¥n(n) = cos mn — =f sin n(z — n)W,(z) ie bw (z)dz. (20) 
0 j 
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The method of successive approximations when applied to this equation will produce 
a converging sequence of solutions. This method is obviously preferable to the direct 
application of the perturbation method, once the equivalence of the results has been 
established, since no minor terms are carried along in the algebraic work, no compli- 
mentary solutions need be added as the integration proceeds,* and the 7; do not ap- 
pear when the function ¢ is evaluated. 

It is of interest to note that when 6;=0 for 7#1, Eq. (20) assumes the form 


2 
b t 
Ys = 05 — f in @ — Oboes, (21) 
0 


and that this equation must generate the elliptic function previously encountered. 
When the method of successive approximations is applied to this equation, the series 
obtained is that one found in the first solution obtained in this paper. This function 
may be obtained more directly by solving Eq. (18) for this particular set of initial 
conditions. 

Perhaps the quickest way to obtain an approximation to the motion for non- 
sinusoidal initial deformation is to be found in the application of a numerical pro- 
cedure using finite differences. Equation (17) lends itself readily to such a treatment 
and the results are considerably easier to interpret than those found by the more 
rigorous integral equation treatment. 

5. The three dimensional problem. If we now allow deflections w normal to the 
plane of u, the procedure of the foregoing sections of this paper leads to the equation 

an? ° 0p fi) 
all 2 2 nS pc 
{1+ flee n + xe leah = (22) 
and to the equation obtained by interchanging g and x in (22). 7 is given by the 
integral on the left side of this equation and x = w/(1+a?r). It follows immediately 
from the similarity of Eq. (17) and that given above that the integral equation 
method previously described will provide the solutions to problems of this nature. 
In particular, however, the motions wherein the string at any instant lies in a single 
plane and wherein each particle describes a quasi-elliptical path is easily determined 
in closed form by considering the deformation expressed in the complex form 


g = ca ¥(ne*™ cos E, 
where w and yu are each real. Equation (22) assumes the form, 
2 ® 2 2 
[1+ 5 f° | ote wirae Fe - (23) 
which, when separated into its real and imaginary parts, implies, 
u'(n) = c/p*(n) 
* When dealing with the differential equations leading to Eq. (8), it was necessary to choose comple- 


mentary solutions to conform to given initial (or other auxiliary) conditions of the problem. In the in- 
tegral equation approach, such conditions are always included in the equations. 
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and 


py +yt 7¥ ~ cy? = 0, (24) 


Here, c is a constant defined by the initial conditions as follows; 
¥(0) = 1, ¥’(0) = 0, u(0) = 0, w’(0) =. 


When c<1+e?/4, these initial conditions lead to a solution of Eq. (24) given by 





, f 1+, — I" ’ 
=]1— (1 —y) sn? ——— €n, a (25 
: seal) : + (on 
” e 
uc ‘ %4s)ds, r= —yP, 
0 4 
where , 
1 es ree 
B = — [VB + €&)? + 326%? + (8 + e*) |, 
Zé” 
1 : 
y = — [V(8 + &)? + 32€%c? — (8 + e*) |. 


¢? 


Note that as c tends to 1+€?/4, y becomes identically unity and the motion of each 


particle is circular. That is, 
ye = ae cos fet 12/4, (26) 


When c>1+e?/4, integration of Eq. (24) yields, 
(8 + 1)(y — 1)2Z? a 
7 (27) 


y? — 1 = ———_—_____—_ 


++8=—@-~ DF 


where 


/1+B yf? —] 
Z=s / —— en, 4/ —— }. 
“Vs 


It is interesting to observe that the string never passes through its rest position for 
values of c different from zero. This follows from the fact that Y never vanishes. 
The function which rigorously defines the transmission of a localized disturbance 
along the string is easily found by considering those solutions of Eqs. (4) which allow 
the function r to assume a constant value. Equations (4) become, under this assump- 





tion, 
0°44’ 024’ Oy’ 0?’ ¥ : 
Pe a ms PP saieee Sms f [(1 + a’r)? — (u’)?}!/2dé = 7, (28) 
o£? On? Og? On? 0 
where 


1i+r 


das a 
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If we now choose 
u’ = f(— — pn), vf = 4(1 + ar)? — (u’)?} 1/2? — 1, (29) 


where 7 is determined by 
1 : d 
JS { [1 + ar]? — | w/(g, 0) |2}/2de = 1, 
us 0 


and where f(£) is non-vanishing in a small region in &, all equations are satisfied. This 
solution is valid until the deformation reaches a fixed point in the string. When this 
occurs, the reflection phenomenon requires a change in 7. This solution is in agreement 
with that found by the linear theory except that p would assume the value unity in 


that theory. 
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A NEW METHOD OF INTEGRATION BY MEANS OF 
ORTHOGONALITY FOCI* 


BY 


A. A. POPOFF 
Automechanical Institute, Moscow, U.S.S.R. 


1. Introduction. This paper contains a new method of integration which is partly 
graphical, partly analytical.! It permits a simple determination of integrals of the 
form [¢;(x)¢:(x)dx, where ¢,(x) is given graphically and ¢(x) is given either graphi- 
cally or analytically. The method requires the construction of certain diagrams, called 
scales, showing the abscissae of the centroids of certain areas associated with ¢,(x), 
and is based on some properties of the so-called orthogonality foci. Finally, the 
method is applied to interpolation, Fourier analysis, and the evaluation of Mohr in- 
tegrals in the theory of structures. 

2. Definite integrals. Let us consider the integral 


1 
T -f i(x)be(x)dx. (2.1) 





If rectangular cartesian coordinates x, y 
are introduced, the functions @¢,(x), 
¢x(x) can be represented by curves, such 
as in Fig. 1. We now consider a distribu- 
tion of mass along the curve y=¢,(x), 









0<x Si, the mass per unit length in the 
x-direction being ¢;(x). The centroid of 
this mass distribution we shall call the 
orthogonality focus. We shall denote it 
by Fx, and its coordinates by px, fi 
(neither p, nor the total mass Q of the 
system depend on ¢;(x)). We have 
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0; -f oi(x)dx. (2.2) 42 | é by a4 
x 
Q; is also the area under the curve j— —j—1—_} +] + 
y =¢,(x). Since T represents the mass "“——————————1— tin 
moment of the mass distribution about Fic. 1 


the x-axis, 


* Russian manuscript received Feb. 24, 1944. The present condensed version was prepared by Dr. 
L. Bers, Brown University, and Professor G. E. Hay, University of Michigan. 
1 This method was announced in the author's note entitled A new method of graphical integration, 


C. R. (Doklady) Acad. Sci. URSS (N.S.) 38, (1943). 


2 This name is justified by the properties discussed in Section 4. 
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1 t 
pr = —f xo, (x)dx, (2.3) fu = T/MQ%. (2.4) 


px is also the abscissa of the centroid of the area under the curve y =¢;(x). 

The following lemmas can be verified easily: 

(a) If (x) is a linear function, its graph is a straight line and F;x lies on this line; 
F;, can thus be found immediately if p;, is known. 

(b) Lf the interval (0, l) is divided into two parts, the orthogonal foci of the two parts 
and of the whole are collinear. 

These two lemmas permit a graphical determination to any desired degree of ac- 
curacy of the point F,, and hence of the integral JT. The procedure is as follows: 

(a) Theinterval (0,7) is divided into 2” equal intervals? (0,/;), (li, /2), - - - , (e™-1,2). 

(b) Operation (a) divides the region under the curve y=@;(x) into 2” regions. 
We find the centroids 4, (r=1, 
2,--+, 2™) of these 2” regions,‘ 
then combine adjacent pairs of 
regions and find the centroids }, 
(7=1, 2,---, 2") of the 27! 
regions so formed, then combine 
adjacent pairs of these 2™~' re- 
gions and find the centroids é, 
(y=1, 2,---, 2*-*) of the 27-* 
regions so formed, and so on. In 
the final stage, we find the cen- 
troid of the entire region under 
the curve y=¢,(x). In the lower 
part of Fig. 1, we see these cen- ee: wae 
troids in a case when m =2. 1 

(c) A diagram, called the scale 






INTEGRAL CURVE 
= > 








as 
ae | 


SCALE OF $4 (0) 




















of (x), is constructed. The mid- vices 

dle part of Fig. 1 shows such a 

scale. It consists of points a,(r=1, 2, - - - , 2”) vertically above 4, and all at the same 
level, points b,(r =1, 2, - - - , 2"-") vertically above }, and all at the same arbitrary 


level slightly below the points a,, and so on. 

(d) Operation (a) divides the curve y=¢@;(x) into 2” parts. We replace each part 
by a segment of a straight line, as shown in the upper part of Fig. 1, and assume that 
the mass is distributed along these segments rather than along the curve. 

(e) We determine the points A,(r=1, 2,- +--+, 2™) of intersection of these line 
segments with the verticals through the points a,. We then determine the points 
B,(r=1, 2,- ++, 2™-") of intersection of the straight lines joining adjacent pairs of 
points A, with verticals through the points },. This process is repeated, until finally 
we arrive at the final point Fy. 

(f) fix, which is the ordinate of Fy, is determined by measurement; the value T 
of the required integral then follows from (2.4). 


3 Unequal intervals could also be used. 
4 It is to be noted that areas corresponding to negative values of ¢,(x) must be considered as corre- 


sponding to negative mass. 
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3. Indefinite integrals. The graphical construction of Section 2 can be applied to 
the indefinite integral /¢;(x)@:(x)dx, in the following manner (see Fig. 2). An in- 
terval (0, /) on the x-axis is taken, and the construction of Section 2 is applied to the 
integral S-bialx)bx(x)dx, where $;,:(x) =¢;(x) in (0, /;) and vanishes elsewhere. This 
yields a point Fj.,, with ordinate fiz. The point (i, fix.) is then plotted. The con- 
struction of Section 2 is then applied to the integral S'bia(x)bu(x)dx, where 
¢:,2(x) =¢;(x) in (0, 2) and vanishes elsewhere. This yields a point Fix,2, and the 
point (Js, fix,2) is plotted. In this way we obtain the sequence of points (l,, fix,r), 


(r=1,2,---, 2™). Since 


1 l 
Sits = — f :,7(x)bx(x)dx, Es. 1) 
sek 0 


the curve passing through these points is approximately the integral curve, except 


for the constant factor 1/{2,. 
Figure 2 shows this construction for the same functions ¢;(x), $%(x) considered 


in Fig. 1, m again having the value 2. 
4. Some properties of orthogonality foci. (a) If the x-axis passes through the 
point F;,, then f;,=0, and by (2.4) 


l 
f o:(x)oi(x)dx = 0, (4.1) 
0 
i.e., &:(x) and ¢;(x) are orthogonal. It is for this reason that Fj; is called the orthogonal- 
ity focus. 
(b) Let us set $;(x) =z(x). Then, from (2.4), 


1 l 
fin = — f [bx(x) |2dx = 2f; (4.2) 
2: J 0 


where f; is the ordinate of the centroid of the region under the curve y=@¢;(x), 
(0x S/). We shall now prove the following theorem. The curve y=h@,(x), where h is 
a constant, has the least mean square deviation from the curve y=;(x) when 


h = fix/frr. (4.3) 
To prove this, we note that the mean square deviation is a minimum when 


d 
dh 


f [oi(x) — hox(x) |*?dx = 0, 
i.e., when 


l 
f [oi(x) — hox(x) |ox(x)dx = 0, 


0 


l l 
h -f si(adon(a)ae / f [ox(x) }2dx = Sir / fur. 


or 
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(c) Let us draw the horizontal line 6 through the point Fy, (Fig. 3), and then rotate 
8 about Fy, through an angle a@ to a new position 8’. A new curve y=¢/ (x) is con- 
structed such that the vertical distance from ’ to points on this curve is equal to the 
vertical distance from 6 to points on the curve y=¢;(x). We shall now prove that 


f $f (x)b(a)dx = f $il2)be(x)de. (4.4) 


We have $/ (x) =¢;(x)+(p,—x) tan a, where, it is recalled, p, is the abscissa of Fy,. 


Thus 


l l 1 
f ob! (x)di(x)dx = f o:(x)b.(x)dx + tan « f (px — x)bx(x)dx. 
0 0 0 


The last integral vanishes, by the definition of px, and the desired result is obtained. 
It is to be noted that the two 
curves have a common point J, about 
which the curve is “rotated.” 
(d) Let us consider the case when |\ 


¢;(x) is linear in each of the intervals iH n - 
(0, 1) (h, 2, so that its graph is a Pjor ||! iv HNN -P 
broken line abcd (Fig. 4). We shall Fix ~ iy NOT 


now show that, if ab is rotated about a 











point M on ab to a new position. a'b’, Doo . p 
then Fy, is unchanged if cd is rotated Ni 
about a certain point N on cd in sucha tH 
way that bc=b'c’. The points M, N 

are called conjugate foci. To prove Pk 





this theorem, we use Fig. 4, in which 











A, Az, Ai, A? are points leading to 1 
the determination of Fy, following 
the procedure laid down in operation Fic. 3. 


(e) of Section 2. From the three pairs 
of similar triangles Mbb’ and MA,Aji, Nec’ and NAsAd, Fy.AiAi and FyA,A}, we 


have 


bb’ By cc’ Bo A,A i ai 











; = 
, 


A 1A 7 Bi- 41 AsAz _— B2 AsA 2 a2 


Since bb’ =cc’, a value for A,;A{ /A2A? can be determined from the first two equations. 
Substitution of this value in the third equation yields 


ayy2(1/B2) + avyi(1/B1) = a1 + ae. (4.5) 


Thus 62 is uniquely determined by (;; hence N is uniquely determined by M. 
It is easily seen that, if rotations of the above type are carried out about conjugate 
foci, and if ¢/ (x) denotes the function the graph of which is a’b’c’d’, then 


i 1 
f él (a)b(a)de = f $i(x)Gu(x)dzx. (4.6) 
0 0 
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(e) Let us denote by Fy, and Fi, 
i the orthogonal foci of the function 
: ¢i(x) represented by the broken line 
\ abcd with respect to two functions 
\ o.(x) and @:(x). We shall now show 
Pe -2 — that, there is a unique pair of con- 
P - ea. — |r oh pe —s jugate foct M (on ab) and N (on cd) 
M - = if i: such that rotations of ab about M and 
ae ie saad pee | cd about N (with bc=b'c’). leave both 
Mm] | \ Fx, and Fj, unchanged. This follows 
>| \ | from the fact that Fy is unchanged 
| * . if 8; and Be (Fig. 4) satisfy (4.5), and 
\ < incisnmimall se tis Os 
Peete \) Shsetest 
r— PB, —--B, ’ a gine - . a 
ee at \ | (4.5). Since both these relations are 
L <a linear in 1/8), 1/B2, they can be solved 
\d" for unique values of 6; and Ps». 
. ' It is easily seen that, if rotations 
pans of the above type are carried out 


about such conjugate foci, and if 
¢/ (x) denotes the function the graph of which is displaced position of abcd, then 


(4.6) holds and also 
l l 
f o! (x)oi(x)dx = f (x) bil x)dx. 


Conjugate foci can be used widely in graphical computations dealing with stati- 
cally indeterminate structures. 

5. Application of orthogonality foci to the interpolation of curves. Let us consider 
the application of orthogonality foci to the following problem. We are given two func- 


tions or curves ¢;(x) and ¢;(x). It is required to find a straight line y=A-+Bx such 


that the integral 
U -f (6; — (A + Bx) |*oidx, (5.1) 


will have the least possible value. We shall now show that U has the least possible 
value when the straight line y=A+Bx passes through the orthogonality foct Fix, Fit, 
where 6; = xox. 

We set 0U/0A =0U/0B =0, to obtain the equations 


i] 


l l ) 
A f od x + B f Xi dx = f didid x, | 
0 0 0 | 
l l l } 
A f xp.dx + BY x*o.dx = f xpipyd x. 
0 0 0 


Let us consider the functions @o(x) =1, ¢1(x) =x, o2(x) =x*. We have Qo =/, 2) = 31’, 
Q, = 4/3. If the orthogonality foci of d, with @o, 1, 62 are denoted by Fio(po, fro), 
Fis(p1, fir), Fia(p2, fiz), respectively, and the orthogonality foci of @; with ¢, and ¢, 
are denoted by Fyx(px, fix), Fis(p1, Fir), respectively, then 
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l l l 
f o.dx = Sirol, f oi. xdx = bfial?, f b.x*dx = $fiol*, 
0 0 0 


£ Gididx = fu fe oidx = firfrol, (5.3) 


ig oioidx = af gi.xdx = Ffiifrl’. 


Thus (5.2) can be written in the form 


2fual 
pe Be as Pee 
2fxo 3fia 








= fit. (5.4) 


Since px, p: are abscissas of orthogonality foci, by (2.3) we have 


l ) 

Pa = f xpd x /f o.dx = — 
0 ko | 

PL = f x*p,.d x / f xo,.d x= 2fual ‘ 
0 0 3 fia } 


Thus (5.5) take the form 


(5.5) 





A + Box = fi; A + Boi = fir, 


whence it follows that A and B must be such that the straight line y= A + Bx passes 
through the orthogonality foci Fy, Fiz. 

In order to construct the straight line y=A-+ Bx which is such that U has the 
least possible value, we can proceed as follows: 

(a) Scales are constructed for do) =1, ¢1=x, ¢2=x*. These are as shown in Fig. 5 
when the interval (0, /) is divided into 8 equal parts (m= 3). 

(b) Ascale is constructed for ¢,. If the function ¢,; is given, this can be done ana- 
lytically. In any event, it can be done graphically using the scales in Fig. 5, since it 
involves integrals of the forms /¢.dx, [oixdx. 
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(c) Ascale is constructed for ¢;=x@,. This also can be done graphically by means 


of the scales in Fig. 5. 

(d) The foci Fy, Fi are found following the procedure outlined in Section 2. The 
straight line through Fy, Fy is the required line. 

6. Graphical harmonic analysis. Orthogonality foci can be used to obtain the 
Fourier series expansion of a function which is given either analytically or graphi- 
cally. If @(x) denotes the function, its expansion into a Fourier series of sines or cosines 


will involve the integrals 


af" _ NUx 
a, = —f (x) sin ; dx (n = 1,2,---), 
0 





(6.1) 





> 
II 


2 : nT Xx 
—f (x) cos dx (n =0,1,---), 
0 


where a,, }, are Fourier coefficients. 
Difficulties are encountered if the method of orthogonal foci is applied directly 
to the integrals in (6.1). These difficulties are avoided if we write (6.1) in the form 


2 l 
a= -f o(x)be1(x)dx, 


i) 


~ [ o(2)ba(a)dx — = fo(ayas oo. oe 





2 7? 2 7! 
— f o@ba(ayax -— = f (aes oe ? or 


> 
ll 


where 





WX nTx 
¢si(x) = sin a Gon(x) = 1+ Sin ; (n = 2,3,+--), 
(6.3) 


na X 





den(x) = COS (n =0,1,---). 
By the use of the scale of ¢o(z)=1 (Fig. 5) and the scales of the functions in (6.3) 
(Fig. 6), the ordinates of the orthogonality foci of ¢(x) with these functions can be 
found graphically by the procedure of Section 2. If we denote these ordinates by 
fo, fens fen, respectively, then (6.2) takes the form 

2 
fe1Qe1, a, = § (Seallen ca SoM) (n = 2, 3, tps ), 


a, = 


(6.4) 


Il 


~| ~| 


b, (fenQen ach fo) (n sas 0, 1, gieate ), 


where Qo, Q,,, Q.n are respectively the areas under the curve ¢o(x) = 1 and the curves 
in (6.3) for the interval (0, 1). Now 


21 1 — cos mr 
Q,=—> 24 =1(14 (mn = 2,3,---), 
Tv nT 


XQ = I, Qe, = l, 
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whence (6.4) becomes 


4 : 1 — cos nr 
a, = — fa, —~™ 2] fou (1 ? Hee) ~ fo] (n = 2,3,---), 
T nT 


by, - Zhen od fo) (n = 0, 1, a >. 


7. Graphical evaluation of Mohr integrals. In the theory of structures, the de- 
termination of deflections in bending often requires the evaluation of so-called Mohr 






| 
| 
| 
| 








> 
T t 














Fic. 7. 


integrals, which have the form T=), M,Mdx, where M, is a function of x given 
graphically and M/ is a linear function of x (Fig. 7). 


From Fig. 7, we see that 
, ee ee - 6 
a= Bae + Be rt (7.1) 
Thus 


_ P Map ' é Mpa 
7 -f ai” —(l — x)dx +f M, % xd x. (7.2) 


By use of the scale of x given in Fig. 5, for both of these integrals the orthogonality 
foci F4z and Fg, can be determined graphically. If Mj, and Mj, denote the ordi- 


nates of these foci, then 


T = Mas(Masl) + Mpa(hMoal), 


or 


T = 1(ManMan + MaaMpa). (7.3) 


* 
M%e2 and M$, are called the focal moments. 
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—NOTES— 


THE PRESSURE DISTRIBUTION ON A BODY IN SHEAR FLOW* 
By M. RICHARDSON (Brooklyn College) 


Problems involving shear flow have been studied recently by Tsien? and Kuo.* 
The purpose of the present note is to point out that the pressure distribution on an 
infinite cylindrical body immersed in a two-dimensional shear flow can be obtained 
by means of integral equations, at least for a sufficiently smooth contour. A direct 
attack on the boundary value problem for the stream function is avoided. The method 
used is essentially that employed by Prager‘ in the case of potential flow. 

If the undisturbed shear flow is given by the velocity field 


v, = U(1 + ky), v, = 0, (1) 


where U and & are constants, then it has constant vorticity equal to —kU. The con- 
tinuity equation implies the existence of a stream function ¥(x, y) which satisfies the 


Poisson equation 


Vy = kU 
with the boundary conditions 
oy oy 
—=-—y9, = 0, —=v, = U(1 + ky) 
Ox oy 


at «, and Y=c, a constant, on the contour C of the cross section of the cylindrical 
body immersed in the usual position in the flow. Let us set Y= Wo+y1, where 


= u(»+5 ‘) 
Yo = y >? 


is the undisturbed stream function and y, is the disturbance stream function. Then 
Vo =kU, and y; is harmonic in the region E exterior to C, with the boundary con- 


dition 


Vie =c- Vos, 
on C, where the parameter s may be the arc length on C measured from any con- 


venient starting point. 

* Received July 6, 1944. 

1 This note was prepared while the author was a fellow in the Program of Advanced Instruction and 
Research in Mechanics at Brown University (Summer 1943). The author is indebted to Prof. W. Prager 
for suggesting the topic and for valuable advice. 

2H. S. Tsien, Symmetrical Joukowsky airfoils in shear flow, Quarterly of Applied Mathematics, 
1, 130-148 (1943). 

3 Y. H. Kuo, On the force and moment acting on a body in shear flow, Quarterly of Applied Mathematics, 
1, 273-275 (1943). 

4W. Prager, Die Druckverteilung an Kérpern in ebener Potentialstrémung, Physikalische Zeitschrift, 


29, 865-869 (1928). 
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By a well-known theorem of potential theory, we have 


1 re) 1 Oy 1 
¥i(P) — yi(~) = —{ (v: —~ log — — —— leg —) as (2) 
ardc On 


On r A 


where » is the exterior normal, P is a point in EZ, and r is the distance between P 
and a variable point whose range will be clear from the context. We now apply 


Green’s theorem 


Ov Ou 
ff (uV*v — vV*u)dA = — f u ;) ds, 
I c on’ ae 


where n’ is the interior normal and Z is the region interior to C, to the functions 
u=—y> and v=log (1/r), obtaining 


1 re) 1 Oo 1 
ev ff log — dA = f vo log —as -f 
I r c On’ r c On’ r 


Using the fact that 0/dn’ = —0/0n and combining this with (2), we obtain 


1 0 1 
¥i(P) — yi(o ) == fv = tog — as - — f tog — as + = ff log—d A. (3) 











On r 


The first integral in (3) vanishes because Y=c on C, and because 


(ee 

— log —ds 

c On . r 

is the angle subtended by C at P, which is zero since P is outside C. In the second 


integral of (3) we may write —dy/dn=v(s) where v(s) is the (tangential) velocity 


along C. Hence we have 


1 
¥i(P) — Wile) = “a v(s) log — ds + —ff log— dA. (4) 
T 
Let us introduce 


1 
V= f v(s) log — ds. 
Cc T 


Then there exist interior and exterior limits 0V;/0n and 0V,./0n such that 
1 (= a) (s) 1 ("4 oV.\ f Ob ) ; 1 -. 
— = Tv(S), = = v(t) — log — dt, 
2 \on On 2 \on an c On . r 


OV, 0 1 
— = — mv(s) + f v(t) — log — dt. (5) 
Cc On r 


On 











so that 


From (4) and (5), we find that the normal derivative of y at the exterior edge of C 
is given by 


OWie 1 kU 
hata te = $09 +5 f 000 tog — ar = ff toe — aa, (6) 


On, 2r On, 
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where the subscript s indicates the point of C at which the quantity is to be evaluated. 


But 
oy rome OYo 2 Wie 




















ns * On, On, On, 
ak —f 1 kU 2 Sf 
= — vs) - — -— ] — aA. 
ling (s) 2x On, vs 


Therefore, the velocity distribution along C, v(s), satisfies the integral equation 


1 0 1 Ovo kU 0 
v(s) + —f v(t) —— log —dt = —2 Sf log — 4A, (7) 
rvJc Ons r 








On, wT ON, 


or, since the last integral may be differentiated under the integral sign, 


v(s) + =f v(t) poe <n ye ff cos (r, we) (8) 


Tet 





where s and ¢ are points of C, (rs, ”.) is the angle between the direction st and the 
exterior normal at s, r is the distance from s to a variable point p of dA, and (r, m,) is 
the angle between the direction sp and the exterior normal at s. This result reduces 
to Prager’s equation (6a), loc. cit.,5 for the special case of uniform flow, that is, when 
k=0. 

The integral equation (7) or (8) for the velocity distribution on the contour C 
may be solved in general by approximative methods. Knowledge of the velocity dis- 
tribution on C is equivalent to knowledge of the pressure distribution on C. 

Example. Suppose C is a circle of radius a with center at 0. In this case, the in- 
tegral equation can be solved explicitly. We have, cos (r.:, 2.)/rs= —1/2a. It is not 


difficult to show that 
fe] 1 
<f log —dA = — za. 
n I r 


Finally, d~o/dn,= U sin 0+ Uka sin? 6 at the point with polar coordinates (a, @). Hence 
(7) or (8) becomes 


a] 


I 
v(s) = —— — 2U. sin 0 — 2Uka sin? 6 + Uka, (9) 


2ar 


where I' = f,0(t)dt is the circulation. 
For the same example, Tsien (loc. cit., equation 18) finds the stream function 


a? k a‘ 
y= U (+ - “) sin 06 + +(r sin? 6 + —— cos 26) |. 
r 2 2r? 


0 
v(s) = — “ = — 2U sin 0 — 2Uka sin? 6 + 3Uka. (10) 


Hence, 


To reconcile this result with (9), we must observe that we can write [=I'9 +I, where 
Ip and T; are the circulations arising from the undisturbed flow and the disturbance 


5 The difference in sign is due to the fact that our (rs, m.) is the angle supplementary to that so de- 
noted by Prager. 
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flow, respectively. Hence Ip = f,vo.dt where vo is the undisturbed velocity field given 
by (1) and the subscript ¢ indicates the tangential component. By Stokes’ theorem, 
Io=Jf,(curl vo)dA=—UkA where A is the area of J. Hence in our example, 
I) = — Ukra’*. If we substitute this for T in (9), assuming, as Tsien does,® that [': =0, 
then our result (9) reduces to (10). 


6 The author is indebted to Dr. Tsien for pointing this out. He had at first mistakenly supposed that 
Tsien’s result was based on the assumption I =0. 


ON PLASTIC BODIES WITH ROTATIONAL SYMMETRY* 
By C. H. W. SEDGEWICK (University of Connecticut) 


Introduction. The rotational symmetry problem in plasticity was discussed by 
H. Hencky! in 1923. In the present paper some new results are obtained. Furthermore, 
the presentation is different from that used by Hencky. 

In the following discussion, 7 and 2 in the cylindrical coordinate system (r, 6, 2) 
will be replaced by a(r, z) and B(r, z) in such a way that a, 8, 6 form a curvilinear, 
orthogonal system. The line element ds will be written in the form 


ds* = A*da*® + B*dB? + r°dé@, 


where A and B are functions of a and 8. Furthermore, if the angle between the curve 
B=const. and the direction of increasing r is denoted by y, we will have 


Or Or 


— = Acosy, — = — Bsiny, (1) 
0a 0B 
Oz : Oz 
—=A siny, — = Bocosy. (2) 
0a 0g 
From these, we get 
OA Oy OB oy 
— a a Be, (3) aneee Si me (4) 
op 0a Oa op 


The stress components will be designated by @aa, 688, 70, Tas, Fa, Ts. In the prob- 
lem under discussion, a9 =@g9 = 0. 

1. Lines of principal stress. Along the lines of principal stress, ¢43=0. In this case 
the equations of equilibrium? reduce to 

* Received December 5, 1944. This paper was written during the summer of 1944 while the author 
was a student in the Program of Advanced Instruction and Research in Mechanics at Brown University. 
The author wishes to express his appreciation to Dr. W. Prager for suggesting the problem and for valu- 
able criticisms. 

1H. Hencky, Uber einige statisch bestimmte Fille des Gleichgewichts in plastischen Kérpern, Zeitschr. 
fiir angew. Math. u. Mech. 3, 241 (1923). 

2 A. E. H. Love, The mathematical theory of elasticity, 4th edition, Cambridge University Press, 1934, 


p. 90. 
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1 oss O o# 9 
— aa) | — — — (In B) - — (1 = (0, 
als ave | rs A aw 











fg ] ©. Da» St deems 
ima a B aB sities 


On simplifying, these become 


OC aa 
— + (Cee — 768) = — * (In B) + (¢aa — on) = —(Inr) = (S) 
Ca 
Ooge 
ia os) = 5g (IBA) + (00 — on) = 5 ann) = 0. (6) 


Assuming the Tresca yield condition, we have ¢aa—@ss = 2k, where k is constant. 
Furthermore, in the so called “fully plastic state,” oo must be equal to either Gaa 
or ogg. Let us assume first that o9 =@aa. Writing Gaa +s +0 = 30, we have Gaa =O 
=o¢+2k/3, o33=0—4k/3. From (5) and (6), we then get 


Oc re) 0a 
— + 2k — (In B) = 0, (7) ~~ 2k = [In (Ar)] = 0. (8) 
Oa Oa 0p 
Elimination of o furnishes 
fs ke 
— |ln (ABr)| = 0. 9 
lin (450) 0) 


If above we had assumed that oo =og3, we would have obtained 


7 op [in(B)] = 0, (79) 8 op nd) =0. — (8*) 
aa da ap ap 
These also lead to Eq. (9), the solution of which is 
A Br = ef (=) ¢0(8), (10) 
Let us define a’ and B’ by? 
da’ = eda, dp’ = e* dp. (i1) 


This transformation merely relabels the families of surfaces a=const. and 
8B =const. 

Now, the volume bounded by the surfaces a, a+da, 8, 8B+d6, 0, 0+d8 is equal 
to A BrdadBd@ = A'B'rda'dB'dé. 

Substituting for da’ and dB’ from (11) and making use of (10), we get 


A'B'r = 1. (12) 


Thus, the volume contained between the co-ordinate surfaces aj , az ; By , Bz ; 61, 02 


is given by 
r oi i da’ dp’ dé = (a? i ay )(B2 _ BY ) (82 _ 6;). 


3 aW. W. Prager, Theory of plasticity, mimeographed lecture notes, Brown University, R. I., 1942. 
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It follows that if the differences a? —a; , 8? —8/ , 62—0, are kept constant for succes- 
sive co-ordinate surfaces, the resulting volumes will be equal. This result is analogous 


to that obtained by Boussinesq‘ in the plane problem. 
Dropping the primes for the sake of simplicity we may construct a solution by 


setting y =g(@). 
From (3), A is then seen to be a function of a alone. We set A =¢’(a), and obtain 


from (4) B=¢g’+h(@). The first equation (1) leads to 
or ¥ , 
r = ¢cos g + I(8), ao — g'¢(sin g) + I’. 


But, according to the second equation (1), 
oe = — ( ) si 
og’ + A) sin g. 
98 g g 


Hence h = —/’/sin g. The condition (12) now takes the form 
l’ 
6 (6¢ = em )o cos g+l/) = 1. 
sin g 


This can be satisfied by setting /=0, g’ cos g=c, @*’ =1/c, where c is a constant. Dis- 
carding constants of integration we thus obtain 


. 3a 
sin g = cB, ¢' =) 
c 
vy = sin (cB), A = 3-2/8¢-18g-2/3, 
B = 31/8¢2/8q1/3[] — ¢2g2]-1/2, r= 3'8e-Miql/3[1 — ¢2g2]1/2, 


Equations (2) now give z =31/8c?/3q/38, Hence 
z/r i cB[1 shes c2p? }-1/2, rp? Lg? = 32/8¢-2/38Q2/8, 


The curves a=const. and 8=const. are thus seen to be concentric circles around 
r=z=0, and radial straight lines, respectively. 

In the above example, it may easily be verified that, corresponding to a set of 
equidistant values of a, 8 and @, the resulting volumes will be equal. 

By substituting the value for B above in (7) and integrating, an expression for @ is 
obtained. 

2. Lines of maximum shearing stress. Along the lines of maximum shearing stress, 
Cap =k and Can =Osg=0. Oo will be equal to either o+k or o—k. Let us assume first 
that oo =o-+Rk. In this case, the equations of equilibrium (2) are 


le a 





|= ve = (4%) . ie 
nha” 6 6h C CUS 8 Uc 





== += = (8%) | + a oe 
imna”:©6 °C CLAD B. a !,lUChe. 


These reduce to 


4 J. Boussinesq, Lois géométrique de la distribution des pressions, dans un solide homogéne et ductile 
soumts @ des déformations planes. Comptes Rendus Ac. Sci. Paris 74, 242 (1872). 
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Oo k or 0A te) 
a [4 — + 2ar | — k—(inr) = 0, 
og 0a 





da ABrL ap 

00 k or 0B o 

— | 3 + 287 —| - k= dn = 0. 

0p ABr da da op 

Making use of Eqs. (1), we obtain 

~ +e sin y- 27] - bane) =0, (13) 
da r 0a 0a 
~ + +f = cosy +22] — k= nr) = 0. (14) 
aB r ap ap 


Eliminating o, we find 


0[—A sin fe] dTBco re) 

ot eee ed ea 7 4 27]. 

0p r da da r 0p 
Carrying out the differentiations and subsfituting for 0A/08, 8B/da, dr/da, Or/dB 
from Eqs. (1), (3) and (4), we obtain 








O*y 1 oy oy AB 
+—[4 cos 7 = — B sin y— | — “cos 2y = 0. (15) 
dadB = 2r 0B 0a 4r? 

We may remark that as r->«, Eq. (15) reduces to that governing the case for 
plane strain, i.e., 0*y/0adB =0. 

It is easily seen that the only solution of equation (15) having two orthogonal 
families of straight lines occurs when y =45°, i.e., when the two families of straight 
lines are inclined at an angle of 45° to the axis of symmetry. This result was obtained 
by Hencky.’ 

If we had assumed above that o~@=oa—k, our equilibrium equations would re- 
duce to 





r 


do —A oy ts] 
—-+ ‘| siny — 2-"| + #= nr) = 0, (13*) 
Oa da 0a 





+ # -* +222] pactiad o 6 (14*) 
po - cos Y Pr Pr nr) = 0, 


which also lead to Eq. (15). 
Let us assume a solution of Eq. (15) in the form y =f(a) +g(8). The equation then 


becomes 
ar[{A cos (f + g)}e’ — {Bsin (f+ g)}f’] — AB cos 2(f + g) = 0. 

Substituting for A and B from relations (1), we get 
dr 8 [cos 2(/f + or 8 [cos 2(f + 
ed ra a a) ka a aie 
da OB 0B da r 

A solution of (16) is given by r=C cos 2(f+g) where C is a constant. Without loss 
in generality we may set f(a)+g(8)=a—8. We then have 


r 
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C cos 2(a — 8B), 


r 
or 


0a 


— 2C sin 2(a — 8B) = — 4C sin (a — B) cos (a — 8B). 


Equations (1) now furnish 
A = — 4 sin (a — 8), B = — 4C cos (a — 8B), 


and Eas. (2) give 
= — 2C(a + B) + C sin 2(a — 8B). 


It follows that the curves a=const. and 8 = const. are cycloids tangent to the lines 


° 
~ 


r=C and r= —C, respectively. 
From Eqs. (14) and (15), we are able to determine o. Substituting our values 


for y, A, B, r and integrating, we find that 
¢ = 4k(a+ 8B) + &l1n [1 — sin 2(a — 6) | + const. 

Another solution is obtained by setting f(a)+g(8)=a—£, as before, and substi- 
tuting r=e*+8d where 6=¢(a—§8) is a function yet to be determined. After making 
these substitutions and carrying out the differentiations, we get 

[¢? — ¢’?] cos 2(a — B) — 2¢¢’ sin 2(a — B) = 0 
which is satisfied by 
¢ = C[cos (a — 8) + sin (a — 8). 
We thus have 
r = Ce*+8[cos (a — B) + sin (a — B)]. 
Using relations (1) and (2), we find that 
z = Ce*t8[sin (a — B) — cos (a — B)]. 

The curves a=const. and 8=const. are logarithmic spirals which intersect the 

straight lines through the origin at an angle of 7/4. This solution corresponds to the 


solution obtained in 1. 
It is interesting to see that these networks of cycloids or logarithmic spirals, known 


in the case of plane strain, are also admissible in the case of rotational symmetry. 


ON THE TREATMENT OF DISCONTINUITIES IN BEAM 
DEFLECTION PROBLEMS* 


By S. TIMOSHENKO (Stanford University) 


In a note on the treatment of discontinuities in beam deflection problems Mr. 
E. Kosko! attributes to R. Macaulay the method whereby the number of constants 
of integration can be always reduced to two, independently of the number of forces. 
This method was, however, originated by A. Clebsch, and is discussed in his book 
“Theorie der Elasticitat Fester Kérper,” 1862, page 389. In Russia it was called the 
Clebsch method and was widely used in textbooks on strength of materials. It was 
also used in German books. See, for example, A. Féppl, Festigkeitslehre, 5th ed. 1914, 


page 124. 


* Received Jan. 14, 1945. 
1 Quarterly of Appl. Math. 2, 271-272 (1944). 
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ON CERTAIN INTEGRALS IN THE THEORY OF 
HEAT CONDUCTIONS* 


By WILLIAM HORENSTEIN (Math. Tables Project, Nat. Bureau of Standards) 


t a? 

o -{ x—*/? exp (- —_—- vx) dx, 
0 x 
t a? 

v -f 21/2 exp (- —_—— ix) dx 
0 x 


frequently arise in the theory of heat conduction. It is the purpose of this note to 
express these integrals in terms of tabulated functions. 
By simple transformations the above integrals may be written in the form 


° b? 1 
o=2f _exp (- on -—)ar, y= m4. J (2) 


vt ? 


The integrals 


(1) 


Let us consider the integral 


x 5? co 
u -f exp (- a*)? — —=)a -f Fay. (3) 


By obvious transformations (3) may be written successively in the forms 


00 c V4 c b\2 
“= f Fd -f Fan = ~* gas — gros f exp[ -at(x+—) |an 
0 0 2a 0 ar 
- b ” b\? 
—— git . — an f \~* exp | - a? ( + =) Ja 
2a a b/ac ar 
Va as b b\? 
—— ¢ ted +4. en f (1 - —) exp| — a7{A+—] |da 
2a b/ac an? ar 
20 b 2 
— ee f exp | - o(r+—) Jar 
b/ac ar 


a ' pa e2ad b 
“= aha cosh 2ab — ba erf (— - ac) 
¢ 


a 2a 


cs) b 2 
— an f exp | - a? (. + -) Ja, (4) 
b/ac ar 


where erf (x) =(2/y mf, exp (—&)dé. 


In an entirely similar manner one finds 


Va b PF ee b\? 
“u = — ¢—* erf (— - ac) + cia f exp | - a? ( ——) |dr. (5) 
2a c b/ac an 


From (4) and (5) one obtains 


lI 


* Received May 2, 1945. 





i 


) 5? 
in A em oom 
exp ( a*n | dx 


Vr Vr b b 
— cosh 2a) + S| eo erf (— = ac) — ¢7% erf (— a a) - 
2a a Cc Cc 


Accordingly, 


) 5? 
¢ 2f _exp (- a*)\? — =) dy 
1/ vt 2 


Jr Vr = a a 
— cosh 2ab + —| e~* erf [| b\/t — —~) — & erf ( bV/t )]. 7 
= cosh 2ab + mm [ er ( \ —) er ( Jt + Ai (7) 


An expression for y is obtained by differentiating (7) with respect to b in accordance 
with the second Eq. (2). ’ 
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